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ABSTRACT 



Two algorithms for generating a non-homogeneous 
Poisscn process with log-quadratic intensity function 
X (t) = exp (ttQ + a^_t + a.2'^^) implemented into 

computer programs and compared for relative speed, 
core storage requirements and fidelity. By simulating 
several cases of non-ho mogeneous Poisscn processes 
with log-quadratic intensity functions it is shewn 
that the Poisson-decomposition and gap statistic 
algorithm, developed by Professor P.A.W. Lewis, Naval 
Postgraduate School, Monterey, California, and G.S. 
Shedler, IBM Research Laboratory, San Jose, 

Calif ernia, substantially reduces computation time 
from that required by an algorithm that uses a 
time-scale transformation of a homogeneous Poisson 
process. The faster algorithm employs a rejection 
technique in conjunction with a method for simulating 
the nen-homogeneous Poisson process with intensity 
function A(t) = exp ( Yq + t) by generation of gap 
statistics. Although additional core storage is 
required by the Lewis and Shedler algorithm, the 
resulting gain in computing efficiency is so 

significant that it outweighs the memory 
consideration. The experience gained from 

implementing the algorithm has led to several 
possibilities which are suggested for improving the 
efficiency of the Poisson-decomposition and gap 
statistic algorithm. 
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IN TRODUCTION 



Many familiar physical and operational processes are 
well described, in whole or part, by examining their "event 
streams" over time. Some such well-known processes are: the 
flow of traffic through an intersection; the arrival of 
persons at seme service facility such as a bank teller's 
window, a service station fuel pump or a grocery store check 
out counter; and, the arrival of telephone calls or radio 
transmissions at some switchboard or other type of 
communications terminal. Analogous processes abound both in 
nature and in the course of our everyday lives. Ey the 
proper definition of an "event" in these situations, the 
process being observed will be characterized by the 
probabilistic nature of the flow, over time, of the events 
of which it is composed. In the above three examples the 
events could be defined respectively as the arrival of a 
vehicle at the intersection, the arrival of a customer at 
the service facility, and the arrival of the telephone call 
or radio transmission at the terminal. The process may then 
be analyzed by examining the interaction of various event 
streams with different intersection configurations, service 
policies and terminal capacities. 

A common method used to perform such analyses is to 
consider the event streams to be homogeneous. This could 
mean that the expected number of events to occur in any two 
or more time intervals of equal length is the same (simple 
homogeneity) , or that the distribution of the number of 
events occurring in any two or more equal time intervals is 
the same (complete homogeneity) . The homogeneous Poisson 
process is often used as a tool in the analysis of such 
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activities. For very simple systems involving event streams 
the use of the homogeneous Poisson process as a model leads 
to tractable analytical results. Most systems of interest, 
however, are not amenable to purely analytical methods, and 
simulation of the processes by digital means becomes 
necessary . 

The assumption of homogeneity in event streams is often 
a very restrictive one. The "rush hour” phenomenon, 
well-known and abundantly cursed by motorists, provides 
cogent evidence that the modeling of event streams is not 
always well served by the homogeneity assumption. The 
intensity of event streams varies over time for many 
physical cr operational processes. In these cases, purely 
analytical methods must be abandoned almost immediately in 
favor of simulation techniques. (The intensity of the event 
stream is defined to be the derivative with respect to t of 
the expected number of events in an interval of length 
(0,t].) 

The ncn-hcmogeneous Poisson process is often employed in 
the analysis of processes that exhibit gross departures from 
the homogeneous event stream criterion. If these processes 
are to be simulated, it is necessary to first describe the 
nature of the inhomogeneity (i.e. how does the intensity of 
the event stream vary over time?) and then to artificially 
generate event streams that behave in accordance with the 
description. 



Of course there 
inhomogeneities that 
However it is int 
streams that display 
types : 



are infinite variat 
can occur or that 
uitively appealing 
varying intensities 



ions in the types of 

can be construed. 

/ 

to consider event 
of the following 



i) increasing continuously over time; 
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ii) decreasing continuously over time; 

iii) increasing and then decreasing continuously over 
time ; 

iv) decreasing then increasing continuously over time. 

A ncn-hcmogeneous Poisson process that has quadratic 
properties can be manipulated to produce the above-mentioned 
effects. The effective simulation of such a process on the 
computer is the subject of this thesis and has been 

motivated by the work of P. A. M. Lewis, Professor, Naval 
Postgraduate School, and G. S. Shedler, IBM Hesearch 

Laboratory . 



There are of course other types of inhomogeneities, such 
as cyclic variations (time of day effect) , but we do not 
consider them here. They have been discussed by COX 
[Hef. 2] and LEWIS [Ref. 9]. In particular LEWIS [Ref. 9] 
describes a process consisting of arrivals at an intensive 
care unit in a hospital. It is shown empirically that, in 
addition to the time of day cycle, long term flucruat icns in 
the intensity function can be adequately described by an 
intensity function whose logarithm is quadratic. 



LEWIS and SHEDLER [Ref. 11] proposed a new 
generating a non-homogeneous Poisson process with 
stream intensity (rate) function that is of 
exponential polynomial form. (The use of e 
polynomials is natural in this context since an 
function is a positive function.) The new method 
have the virtue of increased efficiency over 
conventional time-scale transformation techni 
implemented on a high speed digital computer. 



method for 
an event 
degree-two 
xponential 
int ensity 
appears to 
the more 
que when 



Efficiency in this context is measured in 



terms of 
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computer memory requirements and computational speed. The 
problem of efficiency comparison is recognized by LEWIS and 
SHZDLER in the final pages of their paper [Ref. 11, p. 15]: 
'•There remains the question of efficiency (of the proposed 
algorithm) . . . for generation of a non-homogeneous Poisson 
process with degree-two exponential polynomial rate 
function, relative to generation via time scale 
transformation of a homogeneous Poisson process." 

After seme brief discussion of the requirements of 
implementing the time-scale transformation algorithm the 
report concludes . . . "We therefore would expect the exact 
method of (the proposed algorithm) to be much faster, 
although at the expense of some complexity of programming." 

The objective of this author has been: to implement 
both the algorithm of LEWIS and SHEDLER [Ref. 11 ] and the 
conventional time- scale transformation algorithm on the IBM 
360/67 Computer System in FORTRAN IV language; to define 
reasonable measures of effectiveness for comparing the two 
algorithms; and to determine which algorithm is the more 
efficient in terras of the measures of effectiveness defined. 



Section II discusses the definition of a non-homogeneous 
Poisson process. It also states some special properties of 
Pcisson processes that are used in the development of the 
algorithms investigated in this thesis, and concludes with a 
general discussion of two basic methods of generating a 
non-homogenecus Poisson process. Section III gives a step 
by step description of the two algorithms that were 
implemented into computer programs. The method of 
implementation is the topic of Section IV. Methods of 
comparing xhe algorithms are presented in Section V. 
Section VI presents conclusions drawn from the comparison 
and makes a recommendation for improving the LEWIS and 
SHEDLER [Ref. 11] algorithm. Other recommendations for 
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further study are also listed in Section 71. Appendixes 
provide additional details that would have been awkward to 
include in the main body of the thesis. Computer program 
listings are provided after Appendix E. 
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II. DEFINITION MD PROP ERT IES OF NON- HO KOGENEO US P CISSON 

PROCESSES 



A. GENERAL 



This section will present basic definitions and 
explanations concerning the concept of non-hotnogeneous 
Pcisson processes. References cited may be consulted if a 
mere in-depth understanding is desired. Only the 
fundamental concepts necessary for understanding the 
specific non-homogeneous Poisson process under consideration 
are presented. 



B. DEFINITICN OF A NON-HOMOGENEOUS POISSON PROCESS 

LEWIS and SHEDLER [Ref. 11] define the non-homogeneous 
Poisson process on the real line as follows: 

1. The number of events in any finite set of 
non-overlapping intervals are independent random variables. 

2. Let A (t) be a monotone non- deer easing 

right -continuous function which is bounded on any finite 
interval. Then the number of events in any interval, e.g. 
(O/^q), has a Poisson distribution with parameter 

A(tQ) - A (0) . 

The function A(t) - A (0) is called, among other things, the 
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mean value function since the expected number of events in 
an interval (0,t] is just A (t) - A (0) . Property 1 is the 
independent increment property of the Poisson process; it is 
basic to the idea of a Poisson process. Figure 1 provides a 
graphical representation of the definition above. 

The above definition insures that 

C A (t) - A(0))/{ ACt^) - A(0)} meets the following 

criteria for an arbitrary function F (t) , to be a valid 
distribution function on the interval (0,tQ), c.f. LARSON 



[ Ref . 


7, ch. 3 ]; 






i) 


0 £ F{t) <_ 1 






ii) 


lim F(t) = 0; lim 
t-o t-t„ 


F(t) = 1, 


0 < t < 


iii) 


F(a) _< F(b) for all 


a < b in 


O 

ft 

o 


iv) 


lim F(b + h) = F(b) 


for all b 


in ( 0 



h ^ 0 

where h > 0 . 



Letting F (t) = { A (t) - A (0) } / { A{t^) - A (0) } it follows 

that if A (t) is absolutely continuous in (O/tg] then 
dF(t)/dt = A(t)/{ A (t^) - A(0)} is a valid density function 
on the interval (0,t^]. The function X (t) is called the 
intensity function (or rate function) of the process and 

X(t) = ^ E{nrunber of events in (0,tQ]} 

= ^ {A(t) - A(0)} 
d 

dt • 
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situation: Events occur randomly in time (i.e. along t-axis) 

1) If X = number of events in fixed interval (t^jt^] 

Y = number of events in fixed interval 

Z = number of events in fixed interval 

S = number of events in fixed Interval (0,tQ]; 



Then X, Y and Z must be independent random variables (note: 

S and Z are not independent because of overlap in the intervals). 

2) Given the monotonic increasing right continuous function 
A(t), the number of events in any fixed interval (t.,tj] 
must have the Poisson distribution with parameter A(tj) - A(t.) 
Thus 



X -- Poisson(A(t2) - A(t^ )} 
Y - Po1sson{A(t^) - A(t 2 )} 
Z ~ Poisson{A(t^) - A(t 2 )} 
S - Poisson{A(tQ) - A(0)} 



Figure 1 - DEFINITION OF A N0N-H0.10GENE0US POISSON PROCESS 

(GRAPHICAL REPRESENTATION) 
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C. THE INTEHSITY PDNCTION 



The most intuitively appealing way to think about a 
non-homogeneous Poisson process is to consider the variation 
in the intensity of the event stream over time. The concept 
of an intensity function whose integral meets the criteria 
of the definition in paragragh B above is essential to the 
modeling and simulation of a non-homogeneous Poisson 
process. When one starts with the existence of A (t) , the 
A(t) is often called the integrated intensity (or rate) 
f unction . 



The intensity functuion reveals the instantaneous rate 
of arrivals (in the event stream) as a function of time. 
(This function must be a positive function.) For example, 
if telephone calls arrive at a switchboard at the rate of 5 
per hour at 0900 (9:00 a.m.), increase to a peak rate of 20 
per hour at 1300 (1 :00 p.m.) , then decrease to a rate of 5 

per hour at 1700 (5:00 p.m.), the intensity function could 

lock something like Figure 2a. Then, by plotting the 
integral of the intensity function over the interval of 
interest (i.e. from 0900 to 1700) it is obvious that a 
monotone-increasing, right-continuous and bounded function 
is obtained (see Figure 2b) . If the assumption is made that 
the arrival stream is a Poisson process, i.e. has 
independent increments, then the number of calls received in 
any chosen interval (e.g. 0900-1000, 1230-1315, 1107-1632) 

is distributed as a Poisson random variable with parameter 
equal to the difference between the integral evaluated at 
the right end point of the interval and the integral 
evaluated at the left end point of the interval. These 
values may be read directly from Figure 2b. Specifically, 
the number of calls received in an eight-hour working day is 
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a Ecissoc random variable with parameter 120. 



Although it is unlikely that telephone calls would 

arrive at a switchboard in accordance with the convenient 
parabolic intensity function of Figure 2a, the example 

serves to illustrate two important points. First, it would 
not be realistic to assume that the arrival rate of 

telephone calls at a switchboard would be constant 
throughout a working day. Some function that describes an 
initially increasing and finally decreasing rate of arrivals 
seems more akin to reality. The importance of being able to 
model a non-homo geneous Poisson process is therefore 
established. 

Secondly, it is obvious that the definition of a 

ncn-homogeneous Poisson process is not always used as a 
starting point fcr modeling operational processes. Event 
streams are usually thought of in terms of their underlying 
intensity functions. The idea of an intensity function 
applies to any model for an event stream, and is not 

specific to a Pcisson process. The further step of modeling 
the physical process as a non-homogeneous Poisson process by 
assuming that the process has independent increments is 
taken either on the basis of empirical evidence or physical 
reasoning. Testing for independent increments in a point 
process is discussed in COX and LEWIS [Ref. 3, ch. .6] and 
LEWIS [Ref. 9]. The main physical reason for assuming 
independent increments, and therefore a Poisson process, is 
that the operational process is the superposition of many 
individual event streams. For instance, in a computer 

center egui^ped with several interactive time-sharing 
terminals, the event stream of users at each specific 
teriinal might be assumed to be an arbitrary point process 
with a certain intensity function. The event stream seen by 
the central processing unit is then the sum total (or 

superposition) of the event streams of the individual 
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terminals and, if there are enough terminals, should be 
approximately Poisson. This property of superposition of 
intensity functions is discussed in the next paragraph. 





Figure 2 - EATS AllD INTEGRATED RATS FUNCTIONS FOE TSISPHONE 

CALLS ARRIVING AT A SWITCH30ARD 
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DECOMPOSITION AND SUPERPOSITION OF POISSON PROCESSES 



To ottain a Poisson process with intensity function 
X (t) = (t) + X 2 i^) ^ superpose two Pcisson 

processes, each of intensity X 2 ('*^)* 



We may decompose the intensity function of any 
stream into two or more component event streams. Howe 
we then superpose the component streams, we will recov 
original type process only if we started with a P 
process. For example begin with a renewal process 
intensity v (t) = and let v ^ ^2 ‘ ^ 

processes with intensities and V 2 are superposed 
resulting process is not a renewal process. 



event 
ver if 
er the 
cisson 
with 
enewal 
, the 



This unique pr 
Poisson processes 
intensity functio 
properties of it 
the method used 
Statistic Algorith 



operty may be exploited when simulating 
. It permits the partitioning of the 
n to take advantage of any special 
s component parts. This is the basis for 
in the Poisson-Decomposition and Gap 
m discussed in later sections. 



E. TWO EASIC METHODS OF GENERATING A NON- HOMOGEN EO OS 
PCISSON PROCESS 



1 • Time-Sca le T ransfo rm ati on Algo ri thm 

Consider the non-homogeneous Poisson process with 
intensity function (t) > 0, 0<t<tQ, on the interval 

(0,tQ]. The integral of the intensity function is then 
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A (t) and the number of events in the interval is a Poisson 

random variable with parameter A(t^) - A (0) = (or 

A(tQ) = Pq since A (0) = X(t) dt = 0). Now if T* T* 

are events in a unit homogeneous Poisson process (i.e, a 

Poisson process with constant intensity function 

X ' (t) = X' = 1) then A~^ (T*) , ...,A~^(T*) are events in the 

1 n 

non-homogenecus Poisson process. 



This result gives a procedure for simulating a 

ncn-homogeneous Poisson process, starting with a homogeneous 

Poisson process, which is analogous to the probability 

integral transform method of producing random samples from a 

continuous distribution with distribution function P(x) when 

starting with uniformly distributed random variables. The 

latter method is essentially that if the inverse of F(x) can 

be found, then by generating uniform (0,1) variates u^ , 

. . . , u the values 
n 

= F ^ (u, ) , ... , = F ^(u ) 

i i n n 

comprise a sample of variates from the desired distribution. 



The right-continuous monotone increasing function 
^ (t) describing the ncn-homogeneous Poisson process on the 
interval (0,tQ] can be thought of as a distribution function 
which has been "scaled’* by the factor U q • (Since 
AC^q) = Uq ' ACtQ)/yQ = thus A(t)/^^ is a valid 

distribution function on (0,t^].) 



To implement the time-scale transformation 
one can use the following basic result for ho 
Pcisson processes; given that n events have occu 
homogeneous Poisson process over a fixed interval 
those events are uniformly distributed on the int 
prcof of this property is given in PAHZSN [Ref. 14, 
Therefore, if events are generated as a unit Poisso 
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on an interval of length Mq by first obtaining a 
Poisson (Pq) random variate n, and then letting the order 
statistics from a random sample of uniform (0,Uq) variates 
be the times to events in the unit Poisson process, and the 
times to these events are then transformed by the inverse of 
the integrated intensity function, i.e. A (•) r the effect 
is the same as that obtained from the probability integral 
transform method. Figure 3 illustrates this method 
graphically and also provides a flow chart. Note the 
difference however between this procedure and the 
probability integral transform procedure. In generating a 
unit Poisscn process by this method we need a sample whose 
size is random (and could be zero) i.e. Poisson ( • The 
probability integral transform method simply transforms a 
fixed number cf uniform (0,1) variates into variates from 
seme other distribution. 

In Appendix A it is proved that scaling both the 
distribution function F (x) and the interval (0,1) by the 
same factor does not affect the validity of the probability 
integral transform method. 

The unit homogeneous Poisson process may also be 
generated by adding a sequence of unit exponential variates 
(variates frem an exponentially distributed random variable 
with mean = 1) until the sequence of partial sums of the 
random variables first exceeds Uq (see Appendix D) . This 
accomplishes two things. First it provides an ordered 
sequence cf events from a homogeneous Poisson process on the 
interval (0 ,Uq]. Second, it determines a realization of the 
Poisson randem variable N with parameter A(tQ) = Uq . Note 
that given n, the times to events are uniform order 
statistics. 
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Figure 3 - GRAPHICAL P.3PRSSENT ATION OF TTMZ-SCAL2 

TRANSFORMATION METHOD 
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The difference between these methods for generating 
a homogeneous Poisson process should also be noted. The 
first method requires a Poisson variate and ordering of that 
number of uniform random variables. The second requires 
generaticn of independent exponential variates. The second 
method is probably most basic in that it requires only 
exponentially distributed random variables, and these can be 
obtained from uniform variates by the inverse probability 
integral transform, which is just a logarithm. The method 
is, however, not always the most efficient. 



2 . Conditionin g and Orde r Statistics f a 

Pr oc ess 

This method requires the result of a theorem 
sketched by LEWIS and COX [Ref. 3, ch. 2] and restated here 
for convenience and continuity; it is an extension of the 
result on conditioning in a homogeneous Poisson process 
which results in a conditional uniform distribution of the 
times to events. 

Theo r em Assume that a non-homogeneous Poisson process is 

observed for a fixed time (0,t(^], so that the number of 



events in (0,tQ], 



'to' 



0 

has a Poisson distribution with 



parameter A (t ) - A(0) = p. . Then if T. , . . . ,T denote 

U 0 in 

times-to-events for the process in (0,t„], and if N. = n, 
conditional on having observed n (>0) events in (0,tQ], the 
Tj_'s are distributed as order statistics from a sample with 
distribution function 



= A(t) - A(0) 
A(tQ)- A(0) 



This theorem reduces the problem of simulating a 
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noE-homogeneous Poisson process to that of generating a 
Pcisson number of order statistics from a fixed distribution 
function. That is, given an intensity function over an 
interval (a,b], whose definite integral A (t) = /^x(s)ds is 

Si 

bounded and right-continuous on the interval, an ordered 
sample, ficm a population with distribution function 



F(t) 



A(t) - A(a) 
A(b) - A(a) 



a < t < b (1) 



will yield the desired non-homogeneous Poisson process 

defined by the intensity function X(t) on (a,b]. For 

simplicity the interval will hereafter be assumed to have 
its left end point at zero (a = 0) and its right end point 
at some arbitrary, but fixed point t^, (b = t^) . Using this 
(0,t^3 interval results in no loss of generality, and (1) 

becomes identical with the expression in the theorem. 



Many methods exist for obtaining the necessary order 
statistics. The inverse integral transform explained above, 
decomposition of the density function (see LEWIS Ref. 8), 
or the rejection-acceptance method discussed later in 
Section I7-D, are all possibilities. 

For the family of intensity functions addressed in 
this thesis the n on-hcmogeneous Poisson process is obtained 
by a combination of Poisson decomposition, an algorithm of 
LEWIS and SEEDLEH [Ref. 13] for obtaining a non-homogeneous 
Pcisson process with log-linear intensity function, and the 
conditioning and order statistics theorem given above. The 
procedure involves four basic steps. 



1. The intensity function is decomposed into two 
ccmpcnents; i.e. X(t) = _^(t) + A* (t) . 
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2. A gap statistics algorithm is used to generate 
a ncn-homogeneous Poisson process from one of the 
components, A (t) , which is chosen such that it has 
a special structure (log-linear) . 



3. A rejection-acceptance routine is ussd to 
produce a sample from the remaining component, 
X*(t), i.e. a Poisson number of ordered variates 
are generated. This algorithm is described in 
detail in Section III-C. This sample, when 
ordered, becomes a Poisson process also. 

4. The Poisson processes of the component 
intensity functions are then superposed to produce 
the desired non-homogeneous Poisson process. 



Figures 4, 5 and- 6 (Section III) illustrate the four general 
steps of this procedure. 



Note: Theorem 1 provides a second way to generate 

the unit homogeneous Poisson process required by the 
time-scale transformation algorithm previously described. 
First we may generate a variate from the Poisson random 



variable N. with parameter , say 
Cq u 



N. 



-0 



= n . 



Then 



conditional upon N. = n, n uniform order statistics can be 
generated cn the interval. For reasons explained in 
Appendix D, this second method was used in the computer 
program implementation of the time-scale transformation 
method. 
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F. RATIONALE FOR SELECTION OF DEGREE-TWO EXPONENTIAL 
POLYNOMIAL INTENSITY FUNCTION 



The intensity function identifying the non-homogeneous 
Poisson process investigated in this paper is of the form, 

A (t) = exp(ap, + a, t + a^t^) 



where a^, and a 2 real constants. 



This intensity function was selected for three reasons. 
First, an intensity function must always be positive (or 
zero) if it is to be meaningful. The above function is 



non-negative for all values of a^, 



and a, 2 ‘ Secondly, 



this intensity function, by proper choice of constants, can 
be used to represent the four different types of event 
streams raenxioned previously in Section I. And finally, the 



selection of this intensity function 
statistical procedures; (for details, 
p. 30-34) . 



leads to simple 
see LEWIS Ref. 9, 
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III. COMPETING ALG O RIT HMS 



A. GENERAL 



Assuming an intensity function of the form 



2 

A(t) = exp(aQ + + at ) 

three algorithms for generating the corresponding 
ncn-hofflogeneous Poisson process are discussed. These 
algorithms are based on the two general methods presented in 
Section II-E, and the decomposition and superposition 
property of Poisson processes. 

B. TIME-SCALE TRANSFORMATION ALGORITHM (ALGORITHM A) 



1 . S t e_g O ne 

By definition of a non-homogeneous Poisson process, 
the total number of events observed over a fixed interval 
(0,tg) is itself a Poisson distributed random variable, N^ , 
with parameter Ug = A ( t) dt . The first step of the 

algorithm is to determine the value of the parameter Ug . 
Although an explicit, closed-form expression for the above 
integral cannot, be found, a series representation does 
exist. Except for a constant factor, this series 
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representation assumes the form cf the error function or of 
Daiison's integral. A negative value for the coefficient of 
the second degree term, a , yields the error function form: 




whereas a positive value for results in the Dawson's 

integral form: 




In the above expressions K, t^, and t^ are uniquely 
determined by the coefficients aQ, ( 3.2 points 

of the interval over which the intensity function applies. 
A detailed derivation of above relationships is given in 
Appendix C. 



Evaluation of the error function is a FORTRAN 
supplied procedure and requires only that the proper 
arguments be calculated and provided to the FORTRAN FUNCTION 
ERF or CERF [Ref. 15]. Evaluation of Dawson's integral is 
best accomplished through use of the IMSL (International 
Mathematical and Statistical Libraries, Inc.) FUNCTION MMDAW 
[Ref. 6]. The accuracy of the function values calculated by 
these routines is limited only by the precision 
characteristics of the computer. 



2 . Ste£ Two 



Cnee the parameter of the Poisson random variable 

N is determined, a realization on that random variable is 
0 
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required, (The approach is somewhat backwards since it 
first determines hew many events occurred over the interval. 
It then distributes that fixed number of events over the 
interval in accordance with the non-homogeneous Pcisson 
process described by the intensity function. The importance 
of Theorem 1 now becomes evident since it assures the 
validity of such a procedure.) Generation of Poisson 

variates, especially those with large parameter values, is a 
ccmplex procedure in itself if efficiency in terms of 
computer time and memory requirements is desired. This 
problem is discussed later in Section IV-B. For the 
present, assume that the requisite variate has been 
produced, i.e. N = n. 

^0 

3 • Step T hree 

Given that n events have occurred over the interval 
(0,tQ] we then distribute n events along an interval of 
length in accordance with a homogeneous Poisson process. 

Since events in a homogeneous Pcisson process are uniformly 
distributed over an interval (given that n events have 
occurred), this step merely requires that n uniform (0,1) 
variates be generated, ordered, from lowest to highest, and 
then each multiplied by the factor Uq. The values in this 
n-element vector, (Uj| , u^ , ..., u^) , correspond to the 

points plotted on the vertical axis in Figure 3. 

4 . Step Fo ur 

Each event in the homogeneous Poisson process must 

be transformed by the inverse of the integral of the 

intensity function. Letting A (s) ds = A (t) , the inverse 
, 1 . 0 

A (•) applied to each event in the homogeneous Pcisson 

will produce a corresponding event in the 



process , 
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non^homogeneous Poisscn process, i.e. = t, , 

-1 

A (U 2 ) = etc. The difficulty is that since the 

integral of this specific intensity function cannot usually 

be explicitly expressed, the form of its inverse usually 

eludes any convenient computational formula expression. The 

unique position on (0,tQ] the inverse determines for each 

input value can be found to any degree of accuracy desired, 

by iterative, numerical methods. The Ne w ton-Raphson method 

is easily employed and very efficient in the present 

scenario. Its implementation is explained in Section IV-C. 

Since the function A (t) is strictly monotone increasing, 

the inverse function A~^ (u) applied to an ordered sequence 

of input values results in an ordered sequence of output 

values. Therefore, t^ , t^, ..., t are the times cf events 

12 n 

in the non-hcmogeneous Poisson process and the algorithm is 
complete . 



C. TIME-SCALE TRANSFORMATION ALGORITHM, ALTERNATE 
(ALGORITHM A') 



An alternative approach to the time-scale transformation 
method described above is to generate the required 
hcmcgeneous Poisson process by using the fact that in this 
process the random times between events are independently 
exponentially distributed. Thus one generates unit 
exponential variates until their sum exceeds . The 
partial sums give the times to events and the number of 
partial sums less than or equal to is a Poisson (y^) 
variate. Note that the Poisson variate comes out as a 
by-product in this procedure rather than as a pre-product as 
in Step Two of Algorithm A above. 

Although this method combines Step Two and Step Three of 
Algorithm A into a single procedure, it is not necessarily 
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th€ best method. Because it requires the use of the 
additive method of Poisson variate generation, it becomes 
inefficient for Poisson processes with many events (see 
Appendix D) . For this reason Algorithm A instead of 
Algorithm A’ was used when implementing the time-scale 
transformation method into a computer program. 



D. POISSCN-DECOMPOSITION AND GAP STATISTIC ALGORITHM 
(ALGORITHM B) 



It is recommended that the reader refer to Figures 4, 5 
and 6 for a better understanding of the following steps.) 

1 . St 62 One 

The Poisson- decomposit ion and gap statistic 
algorithm begins with an examination of the coefficients of 
the intensity function. By doing so, the intensity function 
is categorized into one of six possible configurations. 
These six cases are discussed in LSHIS and SHEDLSR 

[Ref. 11]. Examples of each case are illustrated in Figures 
7 through 12 located at the end of this section. 

2 . Ste£ T wo 

a. If the intensity function X (t) is- monotone 
increasing or monotone descreasing over the interval (Cases 
I, II, IV and V; see Figures 7, 8, 10 and 11) the intensity 

function is decomposed into two separate intensity functions 
over the same interval. The resulting intensity functions 
are of the fcrm; 

^(t) = exp(YQ + Y2_t) 

and 
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A*(t) 



= A(t) - _A(t) = expCa^ + a^t + " ®^P(YQ + Y2_t) 

It is clear that _A(t) + A*(t) = A(t). 

b. If either of the two cases not covered by 2a 

abcve occurs, A(t) will be monotone increasing (decreasing) 
on the subinterval (0,b] and monotone decreasing 
(increasing) on the complementary subinterval (b,tQ] (see 
Figures 3 and 6) , where b is a unique point within the 
interval at which A(t) has a maximum (minimum) value. By 
dividing the interval properly into two disjoint, contiguous 
subintervals, each subinterval may be treated as explained 
in 2a. Subsequent steps are applied to each of the two 
intervals separately and the results combined. 



3 . Ste£ T hree 



An efficient algorithm for generating a 

ncn-homogeneous Poisscn process with a log-linear intensity 
function (i.e. A(t) = exp (g^ + g^t) ) is presented by LEWIS 
and SHZELEE [Ref. 13]. This algorithm generates the 

non-homogenecus Poisson process through the use of gap 
statistics. (A comparison of the gap statistic technique 
with the conventional integral transform technique is 

discussed in Appendix C.) By judicious selection cf the 
coefficients cf the log- linear intensity function, most of 
the total area under the original intensity function a (t) 
will be contained under the function A (t) . Therefore most 
of the events in the non- homogeneous Poisson process with 
intensity function A (t) can be accounted for by employing 
the gap statistics algorithm on the intensity function 
A (t) . 
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Step 1 - Categorize intensity function 




Step 2 - Decompose A(t) into ^(t) (log-linear) and 
A*(t) = A(t) - A(t). 




Figure 4 - DIAGRA.^ OF P0ISSCN-DEC0;i?0 SITION A^ID GAP 

STATISTIC ALGORITHM 
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0 1 2 

Step 3 - Generate events T^. from log-linear intensity 
function, ^(t), using gap statistic algorithm. 
Events will be ordered. 
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step 4 - 


Generate events s. 

1 


from intensity 


function 




x*(t) = x(t) - x(t) 


using rejection algorithm 




Events will not be ordered. 
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Step 5a 


- Order events from Step 4, 
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0 1 2 

Step 5b - Merge (superpose) events from Steps 3 and 5a, 
Result is event stream T-]> ... from 

\(t) = Mt) + ^*(t) • 

Figure 5 - DIAGRAM (CONTINUED) 
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Figure 6 - FLOW DIAGRAM OF POISSON-DECOMPOSITIO N AND GAP 

STATISTIC ALGORITHM 
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4 . S te£ Fo u r 



All that remains to be done 
ordered samjle of some given size, on 
from the remaining component of the 
function, i.e. X* (t) 



is to generate an 
the interval (0,tQ], 
original intensity 



The number of events in the interval is a Poisson 



random variable N1 



-0 



with parameter 



u' = /^OX* (t)dt. 

0 



Once a realization on this random variable is obtained, i.e. 

N’ = n'. Theorem 1 may be invoked. Since A*(t)/y' is 
to 0 

more easily evaluated than its indefinite integral, the 

rejection technique (explained in Section IV-D) is used to 

generate the n' required variates. The rejection technique 

is not, in general, always an efficient method for variate 

generation unless great care is taken. Yet in this 

deccmpositicn scenario, it will be used to generate only a 

small percent of the total events required by the original 

intensity function ^(t) . The majority of the events will 

be generated by the efficient gap statistics algorithm. The 

efficiency gains should more than compensate for any 

efficiency losses due to the use of the rejection technique. 



5 . S ten Five 

The events produced by Step 3 will be in order on 
the interval (0,tQ]. The events produced in Step 4 will not 
be in order cn (0,tg]. By ordering the events from Step 4, 
(which are few in number compared to the total number of 
events in the non-homogeneous Poisson process) it is 
possible tc superpose the two ordered event streams. The 
merged event streams produce a new event stream from the 
original intensity function X (t) . 
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IV. 



MiGOIlTHM IMPL5ME NTAT I0H 



A. GENERAL 



This section explains the several specific techniques 
that were applied to impleiaent the two competing algorithms 
(Algorithm A and Algorithm B) into FORTRAN computer 
programs. Detailed discussion of the various subprograms is 
avoided since References 11 and 13 and attached program 
listings provide such information. The Algorithms A and B 
have some subprogram requirements in common while other 
subprograms are unique to one algorithm or the other. This 
section will discuss first those requirements common to both 
algorithms, then those needed only by the time-scale 
transformation algorithm (Algorithm A) and finally those 
unique to the Poisson-decompositio n and gap statistic 
algorithm (Algorithm B) . Hereafter differentiation between 
the algorithms and the FORTRAN computer program 
implementation of the algorithms will not always be made. 
The meaning will be clear from the context. 

B. COMMCN REQUIREMENTS 



• Integ ratio n of a Degree -Two E xpone n tia l P clyn omial 
Q.1§L § fixed Inte rval 

Algorithm A requires that the intensity function 
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A (t) be integrated ever a fixed interval in order to 
determine the value of the parameter of the Poisson random 
variable that governs the number of events that are observed 
to occur in the non-homogeneous Poisson process. Algorithm B 
reguires that the intensity function X*(t) = A(t) - \{t), 

be integrated over a fixed interval for the same reason. 
Since 

b b b 

/ A*(t)dt = / A(t)dt - / ^(t)dt 
a a a 



and A(t) has an explicit expression for its indefinite 
integral, i.e. 



/ A(t)dt = exp(YQ) 
a 



exp (Yj_t) 



Y- 



the problem for both algorithms is reduced to computing the 
value for 

b b 2 

/ A(t)dt = / exp(aQ + a^t + a2t ) dt . 
a a 

SOEEOaTIHE HELP employs IMSL FUNCTION MMDAW or the FCRTHAN 
supplied procedure DERF, as appropriate, to perform this 
calculation. Section III-B and Appendix B discuss how this 
computation could also be made using a convergent series 
representation . 



2 • §®H§£s^pn of a Poiss on Vari a te with a Gi ven 

Parameter 



Both candidate algorithms require at least one 
realization on a Poisson distributed random variable with a 
given parameter. The value assumed by the random variable 



is the number of events observed in a specific interval of 
time over which the occurrence of events is governed by a 
given intensity function (i.e. X (t) for Algorithm A and 
X * (t) for Algorithm B) . The nature of most real event 
streams which lend themselves to analysis using a 
non-homogeneous Poisson process model is that they consist 
of a large total number of events. This will be the case 
either if a dense process is observed over an interval of 
shcrt or moderate length or if a "sparse" process is 
observed over an extended interval (see the arrivals at an 
intensive care unit given in LEWIS [ Ref. 9 ]) . The point is 
that both algorithms must be flexible enough to generate 
ncn-homogeneous Poisscn processes that result in high 
numbers cf events occurring. Since the number of events 
occurring over any interval in such a process is a Poisson 
random variable, it becomes necessary to be able to generate 
a variate from a Poisson distribution with a large 
parameter, i.e. large mean. The two most theoretically 
straightforward algorithms for generating Poisson 
distributed variates (the additive and mult iplioative 
methods) become computationally cumbersome as the parameter 
of the random variable increases. Both the additive and 
multiplicative algorithms and the practical deficiencies of 
each are discussed in Appendix D. 

The technique selected to deliver a Poisson variate 
on demand to both Algorithm A and Algorithem B is the Gamma 
method. It is developed and explained in an unpublished 
book by AHRENS and DIETER [Ref. 1]. A paraphrased account 
of the development of this algorithm is included in 
Appendix C. The main advantage of AHRENS and DIETER'S Gamma 
method is that whereas the computer time required for the 
additive and multiplicative methods is proportional to the 
parameter of the Poisson distribution being sampled, the 
computer time required by the Gamma method is proportional 
to the logarithm of the parameter. The SUBROaTINE PVAR 
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employs the Gamma algorithm to return a Poisson variate when 
given a parameter as an input. 



S torage 



any algorithm which simulates a non-homogeneous 
Poisson process must have some mechanism to provide the user 
with information that completely describes the 
realization (s) of the process simulated. Since the specific 
arrangement of the stream of events over the interval is the 
information of interest to the analyst, the location, or 
time of occurrence, of each single event on the interval 
must be stored. Equivalently, the spacings between events 
would completely describe the realizations on a 
non-homogeneous Poisson process. Thus event spacing 
information rather than event location information could be 
stored. (Programs written for candidate algorithms k and B 
both provide the option for the user to demand either event 
location or event spacing information.) When using either 
algorithm, an array large enough to hold location or spacing 
data for each event generated by the algorithm must be 
created. Since the number of events observed in any 
realization of a non-homogeneous Poisson process is itself a 
random variable, the precise size of the array cannot be 
determined a priori . If the programs are to have any value 
for general application, they must be able to accept 
intensity function parameters which will demand large 
numbers of events when simulated. Using the somewhat 
arbitrary assumption that an event stream with an average of 
4500 events is sufficient for most simulation scenarios, a 
fixed array with a capacity for 5000 events is used in the 
programs implementing both of the algorithms. If the value 
of the intensity function integrated over the interval is as 
high as the maximum limit of 4500 (i.e. the number of events 
in the interval is Poisson distributed with mean = 4500) 
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th€n the array of length 5000 allows the Poisson random 
variable to exceed its mean by 7,45 standard deviations 
before an array overflow is encountered. This highly 
unlikely event is of such rare occurrence (less than one 
chance in a billion) that it may be disregarded. However, 
should it occur, the programs will reinitialize themselves 
and generate a new Poisson process. Also an error indicator 
is incremented. This error indicator (lEH) may be written 
on demand. Its value is the number of times the program was 
forced tc abort generating a Poisson process, reinitialize 
itself and start again. 

A 5000 element array is small enough to avoid its 
being an undue memory requirement burden on most operating 
systems, yet large enough to accommodate most 
non-homogeneous Poisson processes of interest. Its choice, 
though arbitrary, was based upon the above two 
considerations. 

C. SPECIAL EEQUIREMENTS OF THE TIME-SCALE THANSFORMATICN 
AL6CEITHK (ALGORITHM A) 



1. Uniform Variates 



As explained in Section III-B, once the number of 

events observed in the non-homogeneous Poisson process is 

established (i.e. N = n) , it is necessary to construct a 

^ 0 

homogeneous Poisson process consisting of n events over an 



interval of length 



‘0 



units. This is easily done by 



generating n uniformly distributed variates on the interval 
(0,1) and then scaling each by the factor . The LLHANDOM 
computer library package developed by LEWIS and LEARMONTH 
[Ref, 12] is a very efficient source of such variates. Once 
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If 



an array of n uniform (0,1) variates is obtaine 
LLHANDOM, each element of the array must be multipli 
the appropriate scaling factor. It is these res 
numbers which must be acted upon by the inverse o 
integrated intensity function to yield the locat 
events in the non-homogeneous Poisson process on 
original interval (0,tQ]. 



2* Sorting of Event s 

To be of much practical value a simulation r 
for a ncn-hcmogeneous Poisson process should orde 
events in the interval from first to last. In Algori 
this ordering may be done before applying the i 
integrated intensity function to the elements i 
homogeneous Poisson process. The monotonic nature o 
integrated intensity function maintains the relative 
of all elements after they have been transformed b 
inverse function (see Figure 3) . Of course the or 
could be dene after the transformations also, 
implementing Algorithm A, the uniform variates were 
by the factor then ordered, from lowest to hi 

before being transformed by the inverse of the inte 
intensity function. 



Ordering of large arrays of numbers is 
consuming operation on the ' computer. There are 
ordering algorithms of varying degrees of sophisticati 
efficiency. Because ordering is unavoidable when 
Algorithm A, selection of an efficient ordering rout 
important. The ordering algorithm used in 
implementation was the W. R. CHDRCH computer center 1 
routine FXS05T which employs Singleton's version o 
partition exchange sort. (A program listing of PXS 
provided in the computer center's catalogue of 1 
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routines.) The PXSORT routine appears to be the most 
efficient ordering routine readily available for the 
purpose. (It is acknowledged that more efficient routines 
specifically tailored to the problem of ordering uniform 
random variates may possibly have been developed that would 
improve the overall efficiency of the program implementing 
Algorithm A.) 

3. Co mp ut atio n of the Tr ansformed Values 

The essence of Algorithm A is the application of the 
inverse of the integrated intensity function to each event 
in the homogeneous Poisson process on the interval (0 ,uq]. 
As menticned before, the intensity function, X(t), under 
consideration does not yield an explicit expression for the 
integrated intensity function, A(t) = F(t) , that must be 
inverted. (Note: F(t) is a '’scaled'’ distribution function). 
Hence neither does a computationally convenient expression 
for this inverse function exist. Numerical methods must be 
employed to transform the position of each event on the 
interval (0,yQ] in the homogeneous process tc its 
corresponding position on the interval (0,tQ] over which the 
simulated non- homogeneous Poisson process is to be produced. 

The Newton-Ra phson technique was used to accomplish 
the transformation. This technique allows for iterative 
approximations which converge to the true transformed 

values. (i.e. tj_, t 2 r ...f tj^; where t^ = F ^(u^)) The 

iterations continue for each value u^ until an estimate t’ 
for its corresponding t^_ is found that satisfies 

|F(t’) - u^ I < ^ f where ^ is a predetermined tolerance 

-5 

level. (The specific value for £ used was Uq x 10 .) 

In the Newton-Raphso n technique, given a function 
h (X) , the objective is to find the solution, x*, satisfying 
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the equation h (x*) = 0, Letting x^ be the initial estimate 
for X* (based upon some reasonable criteria) the next 
estimate for x*, that is, x 2 / can be calculated from the 
fundamental iteration relationship, 

h(xj^) 

^k+1 " " h' (Xj^) • 

This assumes that the function h(*) and its derivative h'(*) 
can be evaluated at xj^. The process is repeated k times 

until I h (xj^) I < ^ ; or, until 1 Xj^ - I ^ • This 

procedure is nothing more than using the first two terms of 
the Taylcr series expansion of the function h(«) to locate 
the x-intercept of the line tangent to h (x) at the point 
h(Xj^). That intercept value becomes the next approximation 
for the roct of h(x) and the procedure is repeated. A 

graphical explanation of the Newton-Raphson method is 
presented in Appendix E. 

Applying the Newton-Raphson method to the present 
problem requires some special modifications. The problem is 
to find a value t. such that F(t.) = u. , where u. is known, 
i.e. to find t. = F (u.) . Direct application of the 

inverse F (•) is impossible since its functional form 
defies all but the most abstract mathematical expression. 
However if a new function (t) = F(t) - u^ is defined and 
if its root t*, determined (i.e. a t* such that (t*) = 0) 
then F(t=!‘) - u^^ = 0, and F(t*) = u^ . Thus t^ = t* is the 
value desired. 

In order to apply the Newton-Raphson method to the 
function G^(t) it must be possible to calculate both G^ (t) 
and its derivative. Since G. (t) differs from F(t) only by a 
ccnsrant, the problem reduces to that of evaluating F (t) . 
This has already been done, as previously explained in 

Section III-B, and merely requires the assistance of 
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SOEEOOTINZ BELP. 



Now 



G5^(t) = F' (t) , 



so 



derivative cf G^(t) returns the intensity function A(t) 



which is easily evaluated for any 
Ne«ton-Eaphson method may be used. 



t. 



taking 
motion 
Clearly , 



the 



the 



It should be noted here that for each u . there is a 

1 

corresponding function G^(t) cn which the Newton-Baphson 

technique must be used. Since each use of the 

Newton-Baphson method may require several iterations to 

arrive at a value t* which satisfies the tolerance 

criterion, it is readily evident that for large n, 

considerable computational effort is required to obtain all 

the values t^ (i = 1 , ..., n) . The number of iterations 

required to arrive at a suitable approximation for each t^ 

is highly sensitive to the accuracy of the initial 

approximation (designated t ) If the initial approximation 

1 

is close to the actual t^, the Newton-Eaphson method will 
converge very quickly. If the initial approximation is 

poor, convergence could be much slower. The procedure used 
to select initial approximations for each t^ will therefore 
have a profound effect upon the overall efficiency of 
Algorithm A. 

Selection cf these initial approximations is done by 
partitioning the interval (0,tQ] into equal length segments. 
The number of segments is equal to min[10, n/4]. The 
function F(*) is evaluated at the end points of each segment 
and these end points, with their corresponding function 
values are stored in an array. This procedure is performed 
by SDBEODTINE BNCHMK. See Appendix E for a graphical 
representation . 



lo find an initial approximation for the which 
corresponds to any given u^^, the array of function values is 
searched until two adjacent function values which bracket u^ 
are located. These adjacent function values uniquely 
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identify the segment of the interval (0,t^] in which the 
elusive t^ is located. Either end point of the segment 
would serve as a good initial approximation 
Ne wton-Raphson method. However, for the purpose of this 
implementation, the end point which yielded the function 
value closer to Uj_ is used for the initial approximation. 

The decision to divide the interval (0,tQ] into 

min [10, n/4 ] segments was based upon empirical results. Of 

/ 

several proposed segmenting schemes that one which resulted 
in the fewest total number of distribution function 
evaluations over several intensity function/interval 

scenarios was chosen. Higher values than n/4 for sparse 
event streams may produce better results. But n/4 gave good 
results for dense streams and adequate results for sparse 
streams. It appeared to be the best compromise as a 

candidate for general usage. (The number of function 
evaluations of Gj_(t) for each t^ averaged between 2.2 and 
2.7 for this partitioning scheme.) 

Alternative methods would be to use for the initial 
approximation for t one of the following values: 

- t^_i «ti) 

■ \-l 

Neither cf these other methods were attempted so it is 
uncertain whether they would be more or less efficient than 
the method which was used. Efficiency differences among the 
three methods are probably not substantial since each will 
usually yield a good first approximation. 
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D. SPECIAL REQUIREBEHTS OF THE DECOMPOSITION ALGORITHM 
(ALGORITHM B) 



1 . Inte nsity Func^on C ategorization 



The decomposition routine selected depends on which 
of six possible shape categories the intensity function 
X (t) falls into (c.f. Figures 7 through 12) . An 

examination cf the constants cig , and 0.2 in the intensity 
function and the interval (Oftg] over which the 

non-homogeneous Poisson process is to occur will uniquely 
identify the category of the intensity function. A thorough 
discussion of this procedure is presented in Ref. 11, and is 
not reproduced here. Implementation of the category 
identification procedure requires a lengthy sequence of 
decision statements within the computer program. 



2* Selection of the Imbedde d Log-Linear Inte nsity 
Function (s) 



Once the intensity function has been categorized it 
must be decomposed in accordance with a complex scheme 
described in Ref. 11. The objective of the decomposition 
scheme is to fit a log-linear curve (or curves) completely 
underneath X (t) on the interval (i.e. X (t) < X (t) , 0<t<tg) 

in 'such a way as to maximize the area under the log-linear 
curve (s) . This is done by partitioning of the interval 
(0,tg] into subintervals (0,b] and (b,tg] if necessary, and 
then by proper selection of the coefficients Yq and for 

the log-linear function(s). These coefficients are 
functions cf the coefficients and in the original 
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intensity function A (t) and of the interval (0,tQ]. As the 
program advances through the categorizing decision 
statements, the proper coefficients are computed for the 
appropriate leg-linear function(s). 



3 . Ga£ St atistics A lgorithm 



The precursor to LEWIS and SHEDLEB [Ref. 11] was a 
paper by the same authors [Ref. 13] proposing a gap 
statistics algorithm for simulating a non-homogeneous 
Poisson process with a log-linear intensity function 
X (s) = exp (3 q + S]_s) . (Reference 11 reviews this algorithm 
in detail.) As previously noted. Algorithm B divides the 
intensity function \ (t) into a sum of two intensity 
functions, _A (t) and X*(t). The new intensity function 
A(t) is chosen to take on the form of a log- linear function 
so that the cap statistics algorithm may be used to generate 
a stream of events from this portion of the original 
intensity function A (t) . The SUBRODTINE HHPP2 implements 
the LEWIS and SHEDLER gap statistics algorithm for the input 
intensity function ^(t) and returns an appropriate event 
stream to the calling program. Since the use of the gap 
statistics algorithm within Algorithm B is the rationale for 
suspecting it to be a more efficient algorithm than the 
time-scale transformation, it is instructive to examine the 
efficiencies gained in using the gap statistics algorithm 
vice a time-scale transformation algorithm on a leg-linear 
intensity function. This question is addressed in Appendix 
C. It was found that use of the gap statistics algorithm 
resulted in a reduction in computer time of approximately 
50^ from that required by the time-scale transformation 
method. 
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. The Beject ion Routi ne 



The gap statistics algorithm has efficiently 
produced an event stream from a non-homogeneous Poisson 
process defined by the log-linear intensity function 
X (t) = exp (yq + Yi'*') • process to be simulated has 

intensity function X ("t) = «xp (ag + + a 2 ^^)* 

therefore necessary to superpose an event stream from the 
intensity function 

A*(t) = X(t) - A(t) 

2 

= exp(aQ + a^t + a2t ) - exp(Yg + Y^^t) 



onto the event stream obtained from the log-linear intensity 
function. 

Once it is determined how many events are to occur 
over an interval with intensity function A*('h)r i.e. 

= n', it is necessary tc select an ordered sample of n' 
variates from a distribution with density function 



f*(t) 



A* (t) 

\ 

tp, 

/ A*(t)dt 
0 



The value n' will be small compared to the total number of 
events in the original non-homogeneous Poisson process. 
Therefore the need for realizing efficiency in generating 
these n' ordered variates is not usually as crucial to 
overall algorithm performance as is the need for efficiency 
in producing the non-homogeneous Poisson process from the 
log-linear intensity, A (t) , in Step 3. However if the 
technique for obtaining these n* ordered variates is 
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extremely inefficient, much or all of the efficiency gains 
realized from Step 3 above could be lost here. (In fact an 
example of such a loss of previously gained efficiency is 
documented in this thesis; see Section 7I-B.) 

The rejection method is particulary useful for 
generating random variates from populations with continuous 
densities that are bounded and which are concentrated on a 
finite interval; this is the case for f*(t) . The rejection 
method and an algorithm for its implementation are presented 
in LEWIS and SHEDLER [Ref. 11]. k geometric argument will 
suffice to exhibit the principle involved. Consider the 
density function f (x) , for a random variable X on the 
interval (a,b) in Figure 13. The maximum ordinate of this 
density is c. If another function g (x) = c* > c is 
constructed then the density function f (x) is enclosed 
within the rectangle (a,0) , (b,0), (a ,c*) , (b,c*) . If a 

point within this rectangle is selected at random it will 
fall either under the density function or above it. If it 
is under the density curve the abscissa of that point is 
accepted as a variate from the population. If the point 
lies above the density curve that point is rejected and 
another point within the rectangle is selected at random. 
The procedure continues until n' variates have been 
produced. The random points within the rectangle are easily 
produced by generating two independent uniform (0,1) 
variates and scaling them properly resulting in a point 
(x,y), where x = a + u']^(b - a) and y = U 2 *c*. Then if 
f(x) > y, X is accepted as a variate, otherwise it is not. 
After n' variates have been obtained, they are ordered to 
yield an event stream for a process with intensity function 
X * (t) . SDEROOTINE REJECT employs this method in the 
program for Algorithm B. 
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The density f(x) is comoletely enclosed by the rectangle 
(a,0), (b,0), (a,c*), (b,c*). 

Procedure : 

1. Generate 2 ’ndependeT': uniform (0,11 variates u-j 
and u^. 

2. Compute: x = a + bu^ , y = c*U 2 - 

3. Plot (x,y). 

4. If the point (x,y) lies under the density f(x) accept x 
as a variate; otherwise go to Step 1. 

Example : In the graphical example above x^ would be accepted 

as a variate whereas x^ would not be accepted. 
(Note: Ideally, c* = c and a and b correspond 
to the late>"al limits of the density "^(x). This 
minimizes ‘^ejection reqion.) 

Figure 13 - THE REJECTION-ACCEPTANCE METHOD 0? VARIATE 
GENERATION FROM AN ARBITRARY DENSITY 
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since the probability that a point on the abscissa 
will not be accepted as a variate is proportional to the 
area within the rectangle that is not under the density 
function, it is advantageous to make this area as small as 
possible. Therefore it is best to set c* = c and set the 
values a and b equal to the end points of the interval over 
which the density function occurs. This is desirable but 
not necessary. Figure 13 illustrates the more general case 
where the rectangle is larger than necessary. One may be 
required to select a c* > c if c is difficult to determine. 
Seldcm would a and b not coincide with the lateral limits of 
the density however. 

It is obvious from Figure 13 and the preceding 
discussion of the rejection method of variate generation 
that it is the "shape" of the density function which is 
critical to the validity of the method and not the fact that 
it integrates to one. Therefore any scaled density would 
preserve the relative shape of the density function. As 
long as the value c* is at least as great as the maximum 
value of the scaled density, the rejection method will yield 
valid variates. The intensity function A * (t) may be thought 

of as a density scaled by the factor ui = /^®A*(t)dt. 

0 0 

Algorithm B uses the intensity function X* (t) rather than 
the density function A*(‘b)/)jQ in the rejection routine. 
This eliminates a division step for every function 
evaluation required by the method. 

5 • SZ ®ni Stre ams 

The event streams from the intensity functions A (t) 
and A* (t) must be merged to produce an event stream from 
A (t) = A(t) + A * (t) . In the present case there are two 
event streams (one long and one short) which are already 
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ordered. This structu 
than a oiore general 
procedure. To do this 
called. 



re allows for superposition rather 
(and more time consuming) ordering 
merging job SUBROUTINE COLATE is 



E. SUMMARY 



This section has presented a general discussion of the 
more significant and interesting procedures used to 
implement Algorithms A and B efficiently in FORTRAN. The 
reader is referred to LEWIS and SCHEDLEE [Ref. 11] for a 
detailed listing of steps for Algorithm B. Program listings 
may also he consulted. These may be found after Appendix E. 
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V 



• M^HOD OF COMPARISON ^ ALGO RITHM S 



A. GENERAL 



Conclusions drawn from the results of comparin 
efficiencies of two computer programs are only as so 
the measures of effectiveness upon which they are 
These measures of effectiveness are always so 
arbitrary, so their selection should be based 
compelling reasoning. 



g the 
und as 
based, 
mewhat 
upon 



Even after reasonable methods of effectiveness have been 
defined, only a finite number of trials (usually a small 
finite number) for a given set of circumstances may be run. 
Therefore, general statements such as ’’this algorithm is 
better than that algorithm" can seldom be made with complete 
confidence. Yet if rhe incomplete information gathered 
reveals certain trends, generalizations hypothesized from 
such trends may warrant a high degree of confidence. 



This 

compared 



section describes how the two 
and the rationale for choosing the 



algorithms are 
method employed. 
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MEASUBES OB EFFECTIVENESS 



1 . Ccmputat ional Speed 

Perhaps the most popular efficiency measure for 
computer programs is their speed of execution. Speed is a 
natural standard as long as computer time remains a costly, 
scarce resource. 

Algorithms A and B were compared in terms of the 
mean time required to generate event streams from a given 
set of intensity functions. It should be noted that since 
the number of events in each event stream is a random 
variable, the time required to generate one event stream is 
also a random variable. 

Several event streams with different intensity 
functions, each having a different expected number of 
events, are produced. Event streams for each of these 
intensity functions are replicated many times by each 
algorithm. The total execution time required for all 
replications is the "speed" measure of effectiveness. 

The number of replications varies with each 
intensity function. For example, an event stream with an 
expected number of events less than 200 was replicated 100 
times whereas event streams with over 1000 expected events 
were replicated only 30 times. In both cases the total 

number of events produced was of the same order of 
magnitude. A priori it seemed reasonable to assume that 
computer time required to execute each algorithm's program 
would be roughly proportional to the total number of events 
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generated. Therefore in designing the experim 

product of the expected number of events times the 
replications (E[N^^] x # reps) was kept roughly con 
order to keep demands on computer time reasona 
control. (Note: Program execution times were isol 

compilation and linkage times for the purpose 
computational speed comparison.) 



2 . Computer Memory Req uirement s 



A second n 
programs is thei 
were written with 
memory as poss 
complexity. Core 
be accomplished i 
programming method 
weakness as a m 
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3 . F id elity 
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distributed with mean equal to the integral of the intensity 
function over the interval. For a Poisson random variable, 
the variance equals the mean. This test, although useful 
for determining if either algorithm varies grossly from 
expected performance, yields no information about hew the 
algorithms are distributing their simulated events over the 
interval. It is necessary to look at discrete segments of 
the interval and examine the distribution of the number of 
events produced in each segment by our simulated 
non-homogenecus Poisson process. 

From the definition of a non-homogeneous Poisson 
process (see Section II-B) we know that the number of events 
observed over any segment of the interval must be a Poisson 
random variable with parameter equal to the integral of the 
intensity function evaluated over the interval. 3y dividing 
into several segments the interval over which the Pcisson 
process is being simulated, the number of events in each 
segment can be observed. Two hypothesis tests may then be 
made for each segment. 



The first test is the dispersion test [Ref. 3, 

p. 158]. This test seeks to ascertain if the number of 

events observed over each segment is distributed as a 

Poisson random variable. Let n , n , ..., n be k 

X ^ K 

observations on the number N of events occurring in a given 

P_ 

fixed segment p, and let n = (n, + ...+n, )/k. If N is a 

^ p i k p 

Poisson distributed random variable, then the staristic 



d 

P 






2 



n 



P 



has a distribution which is approximated by a chi-squared 
distribution with k - 1 degrees of freedom. If the interval 
has been partitioned into m segments, then by simulating the 
nen-homogeneous Poisson process k times, we can perform m 
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hypothesis tests at a selected significance level a. The 
test statistic dp would then be compared to the (1-a) 
percentile of the chi-squared distribution with k - 1 
degrees of freedom. For a large number of hypothesis tests, 
say 1000, we would expect approximately 1000a rejections of 
the null hypothesis, if in fact the Poisson 

distributed . 



For example consider the non-homogeneous Poisson 
process over the interval (0,100]. If this interval were 
partitioned into 100 unit segments the number of events 
occurring in each segment should be Poisson distributed. If 
51 event streams are simulated over the interval, a value 
for the test statistic d may be computed for each segment, 
i. e . 




Assuming a significance level a= .05 and comparing each d 

2 P 

to the critical value = 67.5048, we would count the 

^ 50 , 0 . 

number of test statistics such that d > 67.5048. If the 

P 

number of events in the segments are indeed Poisson random 
variables we would expect, on the average, that 5 out of the 
100 dp values would exceed the critical value. 

In the above discussion of the dispersion test, no 
mention was made of the parameters of the Poisson random 
variables (assuming that they are Poisson) . This is because 
the dispersion test for Poisson distributions is parameter 
independent. Thus it is necessary to perform a second 
hypothesis test to determine if the mean number of events 
per segment is equal to the difference of the integrated 
intensity function evaluated at the right and left end 
points of the segment respectively. For large k the sample 
mean n is normally distributed. Under the null hypothesis. 
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i.e. that the mean of Nq =/jx(t)dt = (where a. and a. are 

^ a.^ ^ D 

the end points of segment p) n is normally distributed with 

mean yp and variance yp/k* The test statistic is 



c 

P 






and the critical value is the (1-^) percentile of the 
normal distribution, ^a/2* Again, there would be one 
hypothesis test for each segment and if we made 1000 such 
tests we would expect to reject the tull hypothesis 
approximately 1000a times, on the average. 



Continuing the example above, it is now of interest 
to see if the mean number of events on each segment agrees 
with the theoretical value determined by the integrated 
intensity function. This time assuming a significance level 
of a = .10, the critical value for each test statistic Cp 
would be Each Cp > 1.6^5 would be counted and 
if in fact tie random variables have the hypothesized means, 
we would expect (on the average) that approximately 10 test 
statistics cut of 100 would exceed 1.645. 



After the dispersion test and hypothesis 
the means have been performed on event streams 
algorithm, the test results may be compared to see 
algorithm appears to have null hypothesis 
percentages which are not consistent with the a 
the hypothesis tests. 



test for 
from each 
if either 
re j action 
levels of 



(It 

the programs 
determine w 
intensity 
satisfactory 



should be noted that durin 
histograms of the event s 
hether or not they ”fi 
functions. Both program 
results by this subjectiv 



g the devel 
treams were 
t” their 
s seemed t 
e measure.) 



opment of 
plotted to 
respective 
o produce 



63 



r 




c. 



OTHEB CONSIDERATIONS 



1 . In te n sity Function C ategor y 



Algorithm B classifies the input intensity function 
into one of six categories and then determines what action 
to take to generate the appropriate event stream. It might 
be expected that the speed of Algorithm B will be affected 
not only by the total number of events to be generated but 
also by the form of the input intensity function. 
Decomposition of the intensity function into four sections 
will probably require more time than if it must be divided 
into just two sections. 

Algorithm A treats all intensity functions alike. 
Therefore, it is important to compare the two algorithms for 
all six possible categories of intensity functions. It is 
possible that Algorithm B could be faster than Algorithm A 
for one category of intensity function and slower than 
Algorithm A for a different category of intensity function. 
The six intensity functions selected for the comparison 
represent all of the categories defined by Algorithm B, and 
are, in- fact, those intensity functions illustrated in 
Figures 7 through 12 of Section III. These categories are 
numbered frcm I through VI and correspond directly to those 
listed in Ref. 11. 

2 ♦ Ev alua t ion of c* Value in Rej e ctio n R outin e 

The importance of reducing the size of the rectangle 
enclosing the intensity function X* (t) prior to beginning 
the rejection algorithm was mentioned in Section IV-D. If 
possible c* should correspond to the maximum value of X*(t) 
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on the interval (0,t^]. For Cases I and II (Figures 7 and 8) 
this maximum value is easily determined since it occurs at 
one of the end points of the interval. For Case III (Figure 
9) the maximum values of \*{t) and X*(t) both occur at the 
point b where the interval (O^tg] has been partitioned into 
two subintervals (0,b] and (b,t^ ]. For Cases IV, V and VI, 
(Figures 10, 11 and 12) the maximum value of A*(t) is not so 
easily determined. LEWIS and SHEDLEfi [Ref. 11, p. 11] state 
that it is not possible to find thd.s maximum explicitly. 
However an upper bound may be determined from the values Ic 
and k where A(t) < k and \ (t) > k on the interval (0,tQ ] (or 
subintervals (0,b] and (b,tg ] for Case VI). Then by setting 
c* = k - k it is certain that X*(t) < c* on the interval or 
sutintervals of interest. Figure 14 which illustrates the 
rejection routine for the intensity functions for Case V and 
VI used in the analysis, reveals that the c* computed in 
such a manner may not be very efficient as an upper bound 
for the rejection algorithm. One would expect some 
reduction in speed efficiency for Algorithm B when c* 
greatly exceeds the maximum value of A*(t). 



^ Toleran ce Level 



Algorithm A requires selection of a tolerance level 
which determines the stopping criterion for Newton-Raphson 
iterations. The speed of the algorithm could be expected to 
be very sensitive to this tolerance limit. By setting it 
too small the comparison would be prejudiced in favor of 
Algorithm B. By not setting it small enough the fidelity of 
the resulting event stream to the true intensity function 
would be impaired. 

It was determined that a tolerance level of 
E[N^^] X 10 ^ occasionally demanded that the computer 
discriminate beyond its significant digit capability in 
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single precision. This resulted in the failure of the 
Newton-Raphson iterations to converge, although in theory 
they would have done so if enough precision had been 
retained . 



,-6 



A tolerance level of E[N. ] x 10 

'^0 

unnecessarily demanding on Algorithm A even though 



seemed to be 
the 

A stopping rule 
that insured 



computer could discriminate at this level, 
fixing £ = E[N^^3 x 10 ^ was the compromise 
reasonable accuracy of the simulated event stream without 
unfairly burdening Algorithm A. 
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A 



c* = 5 



5 

4 



(Acceptance region accounts for less than 
4% of the total area in the rectangle.) 



max A*(t) = 0.3 



A*(t) 




Figure 14a - Graph of Rejection Routine for the Case V 

Intensity Function Analyzed 




Figure 14b - Graph of Rejection Routines for the Case VI 

Intensity Function Analyzed 

Figure 14 - ILLUSTRATION OF RE JECTION-ACCEFT ANC2 REGIONS 
FOB SAMPLE CASE V AND CASE VI INTENSITY FUNCTIONS 
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VI. 



RESO LTS AND CONCLUSIONS 



A. GENERAL 



Computer programs implementing Algorithms A and S were 
compared in terms of efficiency as explained in Section V. 
This section documents the results of the comparison and 
discusses general observations of the author concerning the 
two competing programs. Also recommendations for further 
study concerning Algorithm B are offered. 



B. MEASURES OF EFFECTIVENESS RESULTS 



1 . Spee d 

Table I convincingly illustrates the superiority of 
Algorithm B over Algorithm A when computational speed is 
used as a measure of effectiveness. The column in Tatle I 
labeled "Time A/Time B" is the ratio of the total time 
required by Algorithm A to produce a specified number of 
replications of the given non-homogeneous Poisson process, 
to the time required by Algorithm 3 to produce the same 
number of replications of the process. For all six 
representative intensity functions the P cisscn-decomposition 
and gap statistic algorithm is at least twice as fast as the 
time-scale transformation. For large event streams the 
superiority of Algorithm B is even more pronounced. 
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Table I - COMPaTATION TIMES 
(IBM System 360/61, FORTRAN 



FOR EVENT STREAMS 
IV (G) compiler) 











Discussion 


Best Relative Performance of B 
over A. Result of c^*" = max X*(t); 
dense event stream; 8% use of 
rejection 


c* = max X*(t); dense event stream; 
20% use of rejection routine 


c* = max X*(t); 14% use of 
rejection routine 


c* > max X*(t); 2nd trial with 
c* - max X*(t) improved ratio; 
18% use of rejection routine 


c* > max A*(t); 2nd trial with 
c* - max X*(t) improved results; 
4% use of rejection routine 


c* > max X*(t); 2nd trial with 
c* - max X*(t) improved results; 
24% use of rejection routine 
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Results obtained for 2nd trial with alternate c* value determined by graphical analysis. 



Ihe most obvious disparity revealed in Table I is 
the difference in speed efficiency (of Algorithm E with 
respect to Algorithm A) between the group of Cases 1 , II and 
III; and the group of Cases IV, V and VI. In the former 
group Algorithm B was 5 to 7.6 times as fast as Algorithm A. 
For the latter group Algorithm B was only 2.2 to 2.5 times 
as fast as Algorithm A. 

This difference between the two groups immediately 
made suspect the value c* used in the rejection routine of 
Algorithm B. For the first group c* is set at the least 
upper bound for x* (*) • second group c* is not, in 
general, a least upper bound for A*(t) but is only an upper 
bound. To determine whether or not this difference in the 
"quality" of the c* values could have such an adverse effect 
on time efficiency for Cases IV, V and VI, these cases were 
simulated a second time with new c* values which were 
empirically evaluated by graphical methods. These modified 
c* values were set close to the least upper bounds of the 
respective A* (t) component of the intensity functions for 
each case. Marked time efficiency gains were observed (see 
Discussion ccluran. Table I) indicating that significant 
efficiency losses can be sustained due to the method of 
setting c* values in Cases IV, V and VI. 

Three factors appear to influence the degree of 
speed efficiency gains of Algorithm B over Algorithm A. 
First is the selection of the c* value as discussed above. 
The second factor is the percentage of total events which 
Algorithm E must generate from the rejection routine. For 
the specific intensity functions analyzed, this percentage 
ranged from 4% in Case V to 24% in Case VI. The fewer 
events required from the rejection algorithm relative to the 
number generated from the gap statistics algorithm, the 
greater the efficiency of Algorithm 3. 
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Finally the total number of events in an event 
stream affects relative efficiencies between Algorithms A 
and B. This is because as the number of events to be sorted 
by Algorithm A grows, the less efficient its sorting routine 
becomes. Sorting requirements for Algorithm B are greatly 
reduced due to the properties of the gap statistic procedure 
used. Therefore as total events increase, so does the 
efficiency advantage of Algorithm B over Algorithm A 
increase . 

2 . Com put er M^or^i Requirements 

Algorithm A requires approximately 72,000 bytes of 
core storage. Algorithm B demands about 92,000 bytes. Both 
programs are costly in terms of storage due to the number of 
library routines used and the need to store up to 5000 event 
times. The difference of 20,000 bytes would not be 
important unless computer capacity were being challenged. 
Memory requirements for both algorithms could be reduced by 
reducing the event stream capacity of the programs. A 
reduction of capacity from a maximum of 5000 events to a 
maximum of 2000 events in each program would result in 
storage requirements of 54,000 bytes and 68,000 bytes for 
Algorithm A and Algorithm B respectively. 

3 . Fidelity 

Table II lists the results of six series of 
hypothesis tests of the types discussed in Section V-B. 
Three series of dispersion tests and three series of 
hypothesis tests for the mean were conducted. Intensity 
functions for Cases III, IV and VI were selected for these 
tests . 
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Fcr €ach of the three cases, a large number (500, 
900 or 1000) of each type of hypothesis test was performed 
for a fixed level of significance a. The number of times 
each null hypothesis was rejected was recorded and the 
statistic a = {(number of re jections) / (total tests)} was 
computed for the series of dispersion tests and the series 
of hypothesis tests for the mean. If the null hypothesis is 
true, then a is an estimate of the significance level a. 

From Table II it appears that for the hypotheses 
tests associated with Case III and Case IV, ct gives a very 
good estimate for a, the true significance level. For Case 
VI some disparities are observed. For Algorithm B the a 
value for the dispersion test (i.e. hypothesis test that 
variates are Poisson distributed) is 0.039 whereas the 
significance level is a = 0.01. Algorithm A yields a much 
better ct value of 0.012. For the hypothesis tests for the 
mean, both Algorithm A and 3 yield a values almost twice 
that of a • 

These results, though interesting, are inconclusive. 
For the dispersion test, the statistic 




is only approximately distributed as a chi-sguared random 

variable. It has been shown by FISHER [Ref. 4] that for 

Poisson random variables with low expectations, hypothesis 

tests based on the chi-sguared distribution may produce 

spurious results. Also, we would expect that if the 

distribution of d is only approximated by the chi-sguared 

p 

distribution, that the relative error between the 

approximate and true distribution is greatest in the tails, 
i.e. at high percentile values. Case VI was tested at an a 
level of 0.01. It is also the intensity function with the 



77 



I 




fewest expected number of events of the three intensity 
functions tested. Partitioning of the interval into 
segments produced some segments which had a low expected 
number of events (see the right tail of the density function 
pictured in Figure 12) . The results for Case VI rather than 
indicating any serious flaws in either of the algorithms 
suggest further research into finding better ways of 
hypothesis testing for sparse event streams and highly 
discriminating levels of significance. 

Eesults of the hypotheses tests on Cases III and IV 
indicate that both algorithms simulate event streams that 
conform to the desired hon-homogeneous Poisscn processes. 
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Table II - RESULTS OF HYPOTHESES TESTS FOR FIDELITY 



{IBM System 360/67, FORTRAN IV (G) compiler) 




Accept 

% 

Reject 

«0 




Accept 

“o 

Reject 

«0 


Hypothesis: Variates Have 

a. 

Mean = / ‘^ A(t) dt = p 


Case VI 
a “ .01 
Trials: 1000 


983 

17 

a = .017 
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Case IV 
a = .10 
Trials: 500 


450 

50 

a = .10 




461 

39 

a =.078 


Case 111 
a = .05 
Trials; 900 


1 1 

1 1 

1 1 

1 I CO 

1 1 to 

CO • C\J I o 

^ 1 to 1 . 

CO 1 i 

1 1 II 

1 1 

1 • < a 

1 1 

1 1 

- L . _ 1 




842 

58 

S = .064 






1 1 




1 1 


Hypothesis: Variates Are 

Poisson Distributed 


Case VI 
a = .01 
Trials; 1000 


1 1 

1 1 

1 1 

1 1 

1 1 CM 

CO 1 CM I O 

CO I 1— 1 • 

0> 1 1 

1 1 II 

1 1 

1 1 < a 

1 1 

1 1 




961 

39 

a = .039 


Case IV 
a = .05 
Trials; 500 


475 

25 

a = .05 




476 

24 

S = .048 


Case III 
a = .10 
Trials: 900 


804 

96 

a = .107 




814 

86 

a = .096 
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c. 



GENERAL OBSERVATIONS 



• Frogrammi nq Eas e 

The programming of Algorithm A is simpler than the 
programming of Algroithra B. Because the time-scale 
transformation technique handles all intensity functions 
alike, several different cases need not be identified. 
Algorithm B must categorize the input intensity function 
into one of six categories and then must treat each category 
in a unique way. A rough indication of this difference in 
the degree of complexity is the number of instruction 
statements required by each program. The program 
implementing Algorithm A required 142 instructions. The 
program for Algorithm B required 364 instructions. 

2. Exact M etho d Ver sus App roxima te Me thod 

Algorithm B employs an exact method in generating 
the non-homogeneous Poisson event stream. It is exact in 
the sense that all event times are calculated directly and 
the precision of those calcularicns are limited only by the 
physical constraints of the computer. 

Algorithm A employs a Newt on-Raphson iterative 
method to approximate event times on the interval. The 
stepping criterion for the technique is limited by machine 
precision considerations. Also the stopping rule does not 
give a firm account of the magnitude of error that can be 
expected in the event times. The epsilon criterion is 
applied to the value of the integrated intensity function 
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and not to the abscissa, (or time interval axis) . For the 
integrated intensity function, F(t), it is conjectured that 
given a function value F (t^) = Uj_, if a t = t‘ can be found 
such that IF(t') ~ F(t^)| < £ that t’ is very close to t 
Although this is a reasonable assumption given the 
well-behaved nature of the integrated intensity function, 
the problem of not knowing how close t' is to t^ may tend to 
reduce one's confidence in the algorithm. 

3 • In it i al i zatio n 



The initialization time required for the programs 
implementing the two algorithms has not yet been mentioned. 
Computation time comparisons did not include compilation or 
linkage times required for the two algorithms. These 
initialization times were significantly different, being 
approximately 16 seconds for Algorithm B and only 8 seconds 
for Algorithm A. Should simulation of only one or two event 
streams be required it is likely that Algorithm A may 
actually require less total time than Algorithm B. The 
difference between the initialization times could probably 
be reduced if the programs were to be rewritten in Assembler 
language . 



D. CONCLCSICN AND RECOaMBND ATIONS 



1 . Conc lusio n 

For general usage, the Poisson-decompositio n and gap 
statistic method of LSMIS and SHSDLER [Ref. 11] is clearly 
superior to the time-scale transformation method in 
generating a non-homogeneous Poisson process with a 
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degree-twc exponential polynomial intensity function. 



The main advantages of the Poisson-decomposition and 
gap statistic algorithm are its speed and its qualification 
as an exact method of variate generation. 

The time-scale transformation method enjoys some 
advantage in core storage requirements and in lower program 
initialization time. These advantages are not sufficient to 
overcome the relative deficiencies of much greater 
computation time and uncertainty about the accuracy of the 
approximating mechanism for determining event times. 

2 • R ecommendatio ns 

The full potential of the Poisson-decomposition and 
gap statistic algorithm can not be realized until it 
includes a better method for selecting an upper limit value 
c* for Cases IV, V and VI intensity functions. Algorithms 
for choosing a c* value that will be close to the maximum 
value for the function X * (t) should be developed and 
incorporated into Algorithm B. A single variable search 
technique (such as a Golden Section search or Bisection 
search) for finding the maximum value of a unimodal function 
may be appropriate. 

The computer program for Algorithm B should be 
rewritten in Assembler language in order to reduce program 
initialization time. The program could then be developed 
into a library routine for general use. 



The question of fidelity of 
ncn-homogenecus Poisson process to the 
process should be investigated further. Of 
is the question of how well the algorithm 



the simulated 
true theoretical 
special interest 
simulates sparse 
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event streams. Methods for conducting the dispersion test 
and the hypothesis test for the mean at very high levels of 
significance (i.e. a = .01, .005) for intervals with a low 
mean number of events are needed. 
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APPENDIX A 



PROOF OF THE VALIDITY OF SCALING THE INVERSE DISTRIEOTION 

FUNCTION 



Propositi on : Let T be a continuous random variable with 

distribution function F (•) , such that 

T 



F^(t) = 0 



P ft) - 

F^(t) 



« < t 1 



F^(t) = 1 



t > t, 



Let Y = such that T € (0,t^) and is a real number. 

Then the random variable Y is uniformly distributed on the 

interval (0 , . 

0 



Proof: maps the line segment [0,t^] onto the interval 

[0,1]. If F (•) is strictly increasing on (0,tQ) then 
(•) exists on the interval (0,1) and maps (0,1) onto 
(0,t^). Now, 



F^iy) = P[Y < y] 






< y] 

P[F (T) < ^)] 

Uq 

P[T < 

F } 

^ ^ uo 

JL 

^0 



0 < y < uq 
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Therefore, 



dFy(y) 

dy 






0 < y < u, 



and Y is uniformly distributed on the interval (0, y.) 
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APPENDIX B 



EBEOR FONCTICN AND DAWSON'S INTEGRAL REPRESENTATION CF THE 
INTEGRATED INTENSITY FUNCTION A (t) 



The standard form for the error function is: 



2 

/? 



t 2 
f e du 
0 



Dawson's integral is usually represented as: 

o t 2 

t / s du . 

0 



Both of these integrals may be calculated to any desired 
degree of accuracy by using the exponential series 
expansion , 



oo n 

X r x: 



nl 



= I 



n=0 



2 2 
u — u 

Thus, s and s may be written as 



2 ® , 2,n <» 2n 

u _ r (u ) ^ r u 

n n! n! 

n=0 n=0 
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and 



2 " , 2,n ® , T >n . 2n, 

,-u _ V (-U ) _ V (-1) (u ) 

n n! n! 

n=0 n=0 



respectively. Integrating yields 



t 2 ® , 2n 

/ du = / I 

0 0 n=0 

” t 2n 

= I f }^du 

n=0 0 



oo ^2n+l 

~ (2n+l) n ! 

n=0 



-2 

Multiplying by t results in 



t „2 



^-2 f U j r t 

t / e du - 2. (2n+l)n 



■» ^2n-l 

I 

n=0 



and the series representation for Dawson's i 
revealed. This series will converge for all t, al 
value of the integral increases very rapidly as t 



ntegral is 
though the 
increases . 
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Using the same argument it is easily shown that the series 

t 

representation for -^q ® is 



(-1) t 



n ^2n+l 



nio 



The error function is obtained by multiplying the integral 
by the constant 2//ir ,i.e.# 



o t 2 _ “ , , , n 2n+l 

^ / e'“ du = -^ I 

/? 0 



✓¥ nio (2n+l)n! 



Although the series expansion may be used to calculate both 
functions, there are very efficient computer library 
routines available for computing Dawson's integral and the 
error function. 

The FORTRAN supplied procedures ERF and DERF [Ref. 15] 
evaluate the above function for input values of t. 



The IMSL (International Mathematical and Statistical 
Libraries, Inc.) FUNCTION MMDAN [Ref. 6] evaluates Dawson's 
integral for input values of t. 

It remains to be shown that the intensity function 

2 

X (t) = exp + a^t + a^t } may be integrated over the 

interval (0,1^] using either ERF (or DERF) or MMDAN. 
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Consider , 



A(t-)-A( 0 ) = / X(t)dt = / exp(a- + a, t + at ) dt 
0 0 ^ 



By completing the square of the exponent the expression 
becomes 

2 

A(t^) = / 0 exp|a„ - | exp | (t + 5^) > dt 



( A (0) = 0, since 0 is the left end point of the interval 
over which \ (t) is defined) . 

For > 0 , 



A (tg) = exp < aQ 



Lettingj 



“l \ 
4 u 2 ( 


0 


exp 




u = 




(t + 


ai 

2a2 


du = 




dt 




K' = 


exp 


“o 


4 a., 



af ) 2 



2a. 



and adjusting the limits of integration, the expression 
becomes: 



TC' ^ K' 

A(t„) = / e du = 



/a 2 a 






b 2 a 2 

/ du - / e’^ du 

0 o 



89 



where , 



and 



a = 






The expressicn above can be rewritten as follows: 



A(t^) 







b 2 
j e du 
0 



2 

a 



-2 

a 




0 



du} 



Inside the brackets are two Dawson integral expressions 
multiplied by constants. Outside the brackets is just one 
more easily determined constant. Using HMDAW twice and 
multiplying by the appropriate constants will give the 
desired integral value. 



For < 0, 



A(tQ) = exp ^ Uq 



al ) ^0 



4a, 



/ 

0 



exp 



■l“2l 5S7’ ^ 



= exp 



a. 



“O ■ 4a. 



0 

/ exp 
0 



“1 2 ) 
+ 2 ^) ] j 



dt 



Substituting as before gives 



A (to) =- 



K 



b 2 
/ e ^ du 



a 

K 



2 



- b 2 

J e du 

/? a 



/¥ 0 



a 2 
/ du, 
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and the error function can be recognized within the 
brackets . 



Two error function calculations, a subtraction and a 

multiplication by a constant are sufficient to determine 

A ( ) • 

0 
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APPENDIX C 



COMPASISCN OF GAP STATISTICS ALGORITHM AND TIME-SCALE 
TRANSFORMATION ALGORITHM FOR SIMULATION OF NON-HO MOG ENEO OS 
POISSON PROCESSES WITH LOG-LINEAR INTENSITY FONCTIONS 



LEWIS and SHEDLER [Ref. 13] present two algorithms 
for simulation of non-homogeneous Poisson processes with 
intensity function X (t) = exp(gQ +6^t) . The first algorithm, 
which uses a time-scale transformation is easy to employ 
since the inverse of the integrated intensity function is 
known explicity. The second algorithm uses gap statistics 
to generate event streams. 

Eoth methods are exact. That is, other than the 
physical precision constraints inherent to the computer, no 
approximations are necessary in generating event times. 
However the gap statistics algorithm obviates the need for 
ordering an array of variates. Also, the gap statistics 
algorithm takes advantage of a fast standard exponential 
variate generator, (Random Number Generator Package 
LLRANDOM, Ref. 12), thereby avoiding the time consuming 
computation of taking logarithms. The time-scale algorithm 
must perform both the ordering of variates and the 
computation cf logarithms. 

Both algorithms were programmed and tested for 
computational speed efficiency. The gap statistics 
algorithm was roughly twice as fast as the time-scale 
transformation algorithm. 
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The SOBROOTINE NHPP2 uses the gap statistics 
algorithm within the Poisson-decomposition and gap statistic 
algorithm (Algorithm B) . It is the presence of SOBRCUTINE 
NHPP2 which markedly reduces the ordering requirements of 
Algorithm B from those necessary in Algorithm A. 
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APPENDIX D 



POISSON VARIATE GENERATION 



A. EACKGHOONC 

Ack nowl edgement; The content of this appendix is a 
paraphrase cf portions of Chapter 11 of an unpublished book 
by J. H. AHRENS and U. DIETER [Ref. 1]. 



The Pcisson distribution has the probability mass 
function 



where is the probability that exactly i events are 
observed to occur in a unit time interval, given that events 
arrive at an average rate \ . 

The intervals between events in a Poisson process of 
rate \ are independent and have an exponential distribution 
with mean i/x, i.e. 



The probability integral transform method can be used to 
obtain a sample of exponential variates from a sample of 



f (x) = Xe 



-Ax 



uniform (0,1) variates, 



since 



X 



1 



X 



2 
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A A rate event stream can then be simulated from the 
exponential variates using the following interpretation: 



1st event occurs at x 

1 

2nd event occurs at x + x 

1 2 

etc . 



This property yields two simple methods for generating a 
realization cn the Poisson random variable. These are the 
additive method and the multiplicative method. 



B. THE ADDITIVE METHOD FOR GENERATING A POISSON VARIATE 



The number k of units arriving in a unit interval in a 
Poisson process with rate A must be found. This will be 
the k for which 



Xt + x_ + • • • + X, <1 
12 k — 



and 



X, + X- + • • • + X, + X, , , >1 

12 k k+1 



or equivalently, for which 



and 



-Jin - Jin U 2 Jin £ A 



-Jin - Zn U 2 -•••- Jin - Jin £ A (1) 



(Note: Since u. < 1 Vi , -In u . > 0) 

1 1 



Thus by using uniform variates, computing logarithms, 
adding and counting, the Poisson variate is obtained. The 
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problem with the additive 
large, say 2000, it requires 
a single Poisson variate, 
fast on a computer, the 
relatively slow. 



method is that if \ is very 
considerable effort to generate 
Although addition is extremely 
computing of logarithms is 



C. THE MULTIPLICATIVE METHOD FOR GENERATING A POI3SON 
VARIATE 



Multiplying (1) by -1 gives: 



5-n (u^ ^ 



and 



Jln(u^‘U2 * 






or equivalently 



and 



Ui-Ur 



U, ‘U_ • ♦ • 

1 2 



• • ♦ • 



u. > e 
k — 



-A 






So by generating uniforms, successive multiplication and 
counting, the Poisson variate k is produced. The 

multiplicative method eliminates the need for computing 
logarithms. However trouble is encountered for large 
values since e ^ becomes very small as A increases. 
Underflow problems occur for A values of approximately 175 
and larger on the IBM 360/57 computer system using single 
precision word lengths. The precision problem may be 

overcome by partitioning the Poisson random variable into 
the sum of two or more Poisson random variables. However 
the average number of multiplications and uniform variates 
required is still proportional to A . More efficient 
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methods exist such as the Gamma method. 



D. THE GAMMA METHOD FOR GENERATING A POISSON VARIATE 



(The following discussion is a condensed account of Ahrens 
and Dieter’s presentation; Ref. 1, pp. 1 1-20/ 21, 22) 

In order to obtain a sample Ic from the Poisson 
distribution of mean y , select a positive integer n 
(typically n is a little smaller than y). Then take a 
sample x from the Gamma distribution of parameter n. 

(1) If x > y return a 

sample k from the binomial 
distribution with parameters 
n-1 / y /x 

(2) If X < y take a sample 

j from the Poisson 

distribution of mean y-x and 
return k n + j. 



The sample x simulates the n-th event (arrival) in a 
Poisson process of rate 1. If x > y then there are n-1 
arrivals in the interval (0,x), and each of these has a 
probability of y/x of being below y [Case (1) ]. If 
X < y then the n simulated arrivals are all before y , 
and the sample j indicates the additional events between x 
and y [Case (2) ]. 

(At this point Ref. 1 includes a formal proof of the 
procedure, which will fce omitted.) 



97 



"•‘"l&iiy - 7''^ 









The explicit algorithm relies on an efficient method for 
sampling from the Gamma distribution. (Such a method has 
been implemented as W. R. CHURCH Computer Center library 
routine GAKA by D. W. Robinson and P. A.W. Lewis and is 
documented in Reference 10.) The constant n is 

arbitrarily chosen as n = [0.875 y ]. This avoids the Case 
(1) almost completely if y is large and a simple method for 
sampling from the binomial distribution becomes 

sufficient: the Bernoulli process is simulated in Steps 
8-10 belcw. The algorithm for sampling from the Poisson 
distribution with mean y - x in Case (2) is the simple 
multiplicative method explained above and is employed in 
Steps 3-5. However, if y - x is still large (larger than 
the "cut-cff point" c) the entire procedure is iterated with 
a new auxiliary sample x from a Gamma distribution of mean 
n* = [ 0.875 (y - x) ] etc. 

Alg orit hm - The Gamma ile^od^ Ahrens a nd Di eter 



1. Initialize k -c 0, w y. 

2. If w > c go to 6 (c = 16 was determined 
empirically by Ahrens and Dieter to be 
reasonable) . 

3. (Start Case (2)) . Set p 1 and calculate 

b ^ e”'^, 

4. Generate u and set p p*u. If p<b deliver 

k . 
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5. 



Increase k k+1 and go to 4 



6. Set n [dw] where d = 7/8 (d value was also 

selected empirically) . Take a sample x from the 

Gamma (n,1) distribution. If x > w go to 8. 

7. Set k = k + n, w w-x and go to 2. 

8. (Start Case (1)) . Set p w/x. 

9. Generate u. If u< p increase k k + 1 . 

10. Set n n- 1 . If n > 1 go to 9. 

11. Deliver k. 

Ahrens and Dieter claim that the computation time for 
the Gamma method grows ultimately like a + B In y . 
Computation time for the additive and multiplicative methods 
grow like u . Therefore for the large y values typically 
demanded in the simulation of non-homogeneous Poisson 
processes, the Gamma method is signally superior in terms of 
computation speed. 
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APPENDIX E 



GRAPHICAL PRESENTATION OF NEWTON-R APHSON METHOD 



/ 

A. GENERAL PRINCIPLES 




Figure E-1 



Consider the function F (t) pictured above in Figure E-1. 
The objective is to find, for a known value Uj_ on the 
vertical axis, a unique corresponding value tj_ such that 
F(tj_) = U;l . If an explicit expression for the inverse of 
F(«) exists, the calculation may be made directly, i.e., 

t = F~l(u ). Otherwise, numerical methods must be used to 
i i 

approximate t , 
i 
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Form the new function G^(t ) = F (t ) 

E-2) . In effect, the horizontal axis 
upward a distance u^. Now if a t* can be 
G^(t*) = 0, then F(t*) = u^ and t* = t^, 



- u^ (see Figure 
has been displaced 
found such that 
the desired value. 



Assume that an initial approximation for t* can be mads, 
say t» . If G. (t? ) and G? (t! ) = g. (t? ) can be comouted 

■‘1 I'l' 1111 

we can write 



9i<4> 



G.(tp 



- 






where 

G^ {tl ) with 

t' gives: 

2 



is the intercept of the line tangent to G^ 
the t-axis. Solving the above expressio 



G.(t') 

t * = t * — . 

^2 ^1 g, (t') 



(•) 

c 



at 

for 



tc 



It is evident from the graph that t^ 
the unknown t*. 



is closer 



than t^ 
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Repeating the procedure (Figure E-3 ) , using t^ as a new 
approximation will yield t^ a still better approximation of 
t*. In general, given an approximation tj^ has been 
obtained, a closer approximation "ktl can be calculated by 
the relationship 









It can be shown that successive approximations will converge 
to the value t* provided Gj^(t) satisfies certain criteria 
[Ref. 5, p. i)49,450 ]. 



Since (t) has the form of a distribution function F (t) 
it satisfies all the conditions for convergence except for 
the case where the interval [tj^ , t^] contains an inflection 
point. If the inflection point(s) can be identified, this 
troublesome situation can easily be avoided. by handling the 
function G^ (•) in discrete sections over the discrere 
intervals, (0,t» ), (t» ,t» ), ..., ) where the t'^ 
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9 



(j = 1, . . . , k) are values such that 

inflecticn points of G over (0,tQ]. 

case under consideration at most one i 
be encountered. It is the maximum 
intensity function.) 



G,-(t". ) are all the 
(For the special 
nflection point will 
or minimum of the 



Successive 
lGi(t]l, ) 1 < £ 
enough to t*. 



approximat ions 
at which time 



are computed until 

t' is assumed to be close 
k 



B. INITIAL AFPEOXIMATIONS 



The problem of obtaining a good initial approximation 
for each was solved in the following manner (see Figure 
E-4 below) : 




Figure E-4 
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The interval (0 , to ] was divided into several segments 
(0/T ], etc. Then the values F(Tj) (j = 1,2,...) 

were computed. Each resulting abscissa and ordinate value 
was stored in an array, i.e.. 



t F(t) 



0 


0 




F ( t ^) 


^2 


F(t2) 


• 


• 


^0 


F(to) 



For each value u the array vias searched until two function 
values were located such that F(Xj) < u < • Either 

T. or was then selected as the initial approximation 

for the value t- (i.e. t- ,) such that F(t-) = u • . The 

X. X. ^ X J- X 

function G^{t ) = F(t ) - u^^ was then formed and the 

Newton-Saphscn iterations performed until the epsilon 
criterion was satisfied. 
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SUBROUTINE MTOINV 
PURPOSE 

SIMULATES A NON-HOMOGEN EOUS POISSON PROCESS WITH 
QUADRATIC EXPONENTIAL INTENSITY FUNCTION OVER A 
GIVEN INTERVAL USING A TIME-SCALE TRANSFORMATION 
OF A HOMOGENEOUS POISSON PROCESS. 

USAGE 

CALL MTOINV ( A, A1 , A2 , ELt ER, I S , 1 1 , M» I ER ) 

DESCRIPTION OF PARAMETERS 

A - CONSTANT IN INTENSITY FUNCTION. 

Al - 1ST DEGREE COEFF IN INTENSITY FUNCTION. 

A2 - 2ND DEGREE COEFF IN INTENSITY FUNCTION. 

EL - LEFT END POINT OF INTERVAL. 

ER - RIGHT END POINT OF INTERVAL. 

IS - RANDOM NUMBER SEED. ANY INTEGER WITH NINE 
OR LESS DIGITS. 

II - 0 FOR TIMES OF EVENTS. 

I FOR TIMES BETWEEN EVENTS. 

M - RETURNED VALUE. EQUALS NUMBER OF EVENTS IN 
NON-HOMOGEN EOUS POISSON PROCESS. 
lER - ERROR CODE. MAY BE WRITTEN ON DEMAND. I ER 
EQUALS THE NUMBER OF TIMES PROGRAM WAS 
FORCED TO ABORT SIMULATION AND START AGAIN 
DUE TO EXCESSIVE EVENTS (MORE THAN 5000). 

COMMENTS 

calling program must have a common REGION, BLOKl, 
OF DIMENSION (5000). 

EXAMPLE: DIMENSION T(5000) 

COMMON/ BLOKl/T 

CALLING PROGRAM MUST USE THE FOLLOWING JCL CAROS 

// EXEC FORTCLG, IMSL=DP 
//FORT.SYSIN OD * 

TIMES TO EVENTS OR TIMES BETWEEN EVENTS WILL BE 
STORED IN CELLS T(l) THROUGH TIM). 



SUBROUTINE MTDINV ( A , Al , A2 , EL , ER , IS , 1 1 , M , I E R ) 
DIMENSION TAU(5000), 9RKT(5000,2) 

COMMON /BLOKl/ TAU/BL0K2/BRKT 
CALL OVFLOW 
lER = - 1 
1 lER = IER+1 
BRKTI 1,1) =0. 

BRKTI 1, 2) = 0. 

INTEGRATE INTENSITY FUNCTION 

CALL HELP (A,A1,A2,EL,ER,PAR) 

PARI = PAR 

GENERATE POISSON VARIATE M 

CALL PVAR ( IS, PARI, M) 

IF (M.GT.5000) GO TO I 

SET EPSILON VALUE FOR N EWTON-RA PHSON STOPPING RULE 

EPS = FLOAT (M+1 ) *0. OOOOl 
C 
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GENERATE M UNIFORM VARIATES; ORDER AND SCALE 

CALL RANDOM ( IStT AU,M) 

CALL PXSORT (TAU»ltM) 

00 2 1=1, M 

TAU( I ) = TAU( I )#PAR 

2 CCNTINUE 

DETERMINE INTERVAL PARTIONING SCHEME 
FOR SUBROUTINE BNCHMK 

12 = M/4 

OIV = AMAX0( 10, 12 ) 

IDIV = INT(0IV+2) 

DELTA = (ER-ED/DIV 
ENOL = EL 

COMPUTE ARRAY OF ABSCISSAS AND ORDINATES FOR 
INTEGRATED INTENSITY FUNCTION 

CALL BNCHMK ( A, Al , A2 , ENOL , 0 ELT A , I DI V ) 

BEGIN NEWTON-RAPHSON ITERATIONS 

K = 0 
J = 1 

3 K = K+1 

4 IF (TAU(K).LT.BRKT( J ,2) ) GO TO 5 
J = J + 1 

GO TO 4 

5 IF ( ABS(BRKT( J, 2)-TAU(KJ) .LT.A8S(BRKT( J-1 ,2 )-TAU( K) ) ) 
*GC TO 6 

T = BRKT(J-l,l) 

ADO = ( BRKT(J-1,2)-TAU(K) )/EVAL(A,Al,A2,T) 

GO TO 7 

6 T = BRKT( J, I ) 

ADD = (BRKT(J,2)-TAU(K) )/EVAL(A,Ai,A2,T) 

7 IF (ABSCAOOl.LT.EPS) GO TO 8 
T = T-AOO 

CALL HELP ( A, A1 , A2, EL, T,X ) 

ADO = ( X-TAU(K) )/EVAL(A,Al , A2,T) 

GO TO 7 

8 TAU(K) = T 

IF (K.EQ.M) GO TO 9 
GO TO 3 

9 IF (lI.EQ.O) RETURN 

TIMES BETWEEN EVENTS ARE REQU EST ED-CALCU LA TE THEM 

S = T AU ( I ) 

TAU( 1 ) = TAU( 1 ) -EL 

00 10 1=2, M 
SI = TAU( I ) 

TAU( I ) = TAU( I)-S 
S = S 1 
10 CONTINUE 

RETURN 

END 



SUBROUTINE HELP CALCULATES THE VALUE OF THE 
INTEGRATED INTENSITY FUNCTION OVER ANY INTERVAL (0,T) 
OVER WHICH THE NHPP OCCURS 

SUBROUTINE HELP ( A , A 1 , A2 , EL , ER , XX ) 

DOUBLE PRECISION MMDAW,BB,AA 
Z = SQRT( ABS1A2) ) 

Y = (A1=!'Z) /I 2.*A2 ) 
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AA = Z«EL+Y 
B8 = Z=!‘eR+Y 
CC = A-Al-A 1/ ( 4.*A2) 

CC = EXP(CC)/Z 
IF (A2.LT.0.) GO TO 1 
Q1 = DEXP( AA**2 ) *MMOAW( AA ) 
Q2 = DEXP(B8=«'=«'2)-MM0AW(BB) 
XX = CC«(Q2-Q1) 

RETURN 

1 CC = CC*. 8862269 

XX = CC*(DERF(8B )-DERF( AA) ) 

RETURN 

END 



SUBROUTINE PVAR GENERATES A POISSON (THETA) 

VARIATE, M, USING THE GAMMA METHOD 

SUBROUTINE PVAR (IS,THETArM) 

DIMENSION T(5000) 

COMMON /8L0K1/ T 
K = 0 
C = 16. 

D = .375 

1 IF (THETA. LT.C) GO TO 2 
GO TO 5 

2 U = 1. 

CTN = EXP(-THETA) 

MMAX = THETA+6.*SQRT(THETA) 

3 1 = 1 

CALL SRAND (IS, T, MMAX) 

4 U = U*T( I ) 

IF (U.LT.CTN) GO TO 8 
I = I + l 
K = K+1 

IF ( I.GT.MMAX) GO TO 3 
GO TO 4 

5 NP = INT(0*THETA) 

AN = FLOAT (NP) 

CALL GAMA ( AN, IS,G,1 ) 

IF (G.GT. THETA) GO TO 6 
K = K+NP 
THETA = THETA-G 
GO TO 1 

6 U = THETA/G 
NP = NP-l 

CALL SRAND (IS,T,NP) 

DO 7 1=1, NP 

IF (T(I ).LT.U) K = K+1 

7 CONTINUE 

THE VALUE M IS ASSUMED BY THE POISSON (THETA) VARIATE 

8 M = K 
RETURN 
END 



SUBROUTINE BNCHMK GENERATES ARRAY OF ABSCISSAS 
AND ORDINATES OF THE INTEGRATED RATE FUNCTION. 
THIS ARRAY IS USED TO DETERMINE INITIAL ESTIMATES 
FOR NEWTON-RAPHSON ITERATIONS 

SUBROUTINE BNCHMK ( A , A1 ,A2 , B I NC , CELT A, I D IV ) 

DOUBLE PRECISION AA,36,MMDAW 
DIMENSION BRKT(5000,2) 

COMMON /BLOK2/ 8RKT 
Z = SQRT ( ABS( A2 ) ) 

Y = (Al*Z)/(2.*A2) 
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AA = Y 

CC = A-A1^A1/(4.*A2) 

CC = EXP(CC )/Z 

IF (A2.LT.0.) GO TO 2 

Q1 = OEXP( AA**2)=«'MM0AW(AA) 

C 

DC 1 I = 2,IDIV 
BINC = BINC+OELTA 
BB = Z=<'BINC+Y 

Q2 = 0EXP( B8=J==«^2 ) ^^'MMOAWl BB) 
XX = CC-(Q2-Q1) 

8RKT( 1,1) = BINC 
BRKT( 1,2) = XX 

1 CONTINUE 
C 

GO TO 3 

2 CC = CC*. 8862269 
C 

DO 3 I=2,IDIV 
BINC = BINC+OELTA 
BB = Z*BINC+Y 

XX = CC*(OERF(BB)-OERF( AA) ) 
BRKTU,1 ) = BINC 
BRKTI 1,2) = XX 

3 CONTINUE 
C 

RETURN 
END 



FUNCTION EVAL COMPUTES THE VALUE OF THE 
INTENSI TY FUNCTION 

FUNCTION EVAL (A,A1,A2,T) 

EVAL = A+A1*T+A2*T**2 
EVAL = EXP( EVAL ) 

RETURN 

END 
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SUBROUTINE DEGTWO 
PURPOSE 

SIMULATES A NGN-HOMOGENEOUS POISSON PROCESS WITH 
QUADRATIC EXPONENTIAL INTENSITY FUNCTION OVER A 
given interval using the POISSON-DECOMPOSITION 
AND GAP STATISTIC ALGORITHM. 

USAGE 

CALL DEGTWO ( I S t A , A1 , A2 , E L , ER , 1 1 , Nt I ER ) 

DESCRIPTION OF PARAMETERS 

IS - RANDOM NUMBER SEED. ANY INTEGER WITH NINE 
OR LESS DIGITS. 

A - CONSTANT IN INTENSITY FUNCTION. 

Ai - 1ST DEGREE COEFF IN INTENSITY FUNCTION. 

A2 - 2ND DEGREE COEFF IN INTENSITY FUNCTION. 

EL - LEFT END POINT OF INTERVAL. 

ER - RIGHT END POINT OF INTERVAL. 

II - 0 FOR TIMES OF EVENTS. 

1 FOR TIMES BETWEEN EVENTS. 

N - A VECTOR OF LENGTH 5. N( 1 ) THROUGH N(4) 
CONTAIN NUMBERS OF EVENTS FROM VARIOUS 
COMPONENTS OF THE DECOMPOSED INTENSITY 
FUNCTION. N(5) CONTAINS THE TOTAL NUMBER 
OF EVENTS IN THE NON-HOMOGENE OUS POISSON 
PROCESS . 

COMMENTS 

CALLING PROGRAM MUST HAVE A COMMON REGION, HOLD, 
OF DIMENSION (5000), AND AN INTEGER ARRAY OF 
DIMENSION (5). 

EXAMPLE: DIMENSION T(5000), N(5) 

COMMON/HOLO/T 

CALLING PROGRAM MUST CONTAIN THE FOLLOWING 
ASS IGNMENT STATEMENT: 

M=N( 5) 

CALLING PROGRAM MUST USE THE FOLLOWING JCL CAROS 

// EXEC FORTCLG, IMSL=OP 
//FORT.SYSIN DD * 

TIMES TO EVENTS OR TIMES BETWEEN EVENTS WILL BE 
STORED IN CELLS T(l) THROUGH T(M). 



SUBROUTINE DEGTWO ( I S , A , Al , A2 , E L , ER, I I , N , I ER) 
DIMENSION TIMES(5000), T(5000), N(5), P(5) 
COMMON /MIKE/ TIMES/HOLO/T 
CALL OVFLOW 

INITIALIZE variables 

P(l) = A 
P(2) = AI 
P(3) = A2 
P(4) = 0. 

P(5) = 0. 

C 

DO 1 1=1,5 
1 N( I ) = 0 
C 

C IF RATE FUNCTION IS LESS THAN DEGREE TWO, 
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USE NHPP2 ROUTINE ONLY 

IF (A2.EQ.0.) GO TO 2 
GO TO ^ 

2 CALL NHPP2 ( I S , EL , ER , A, Al , I I , N I , I ER ) 

N(5) = N1 

IF (Nl.EQ.O) RETURN 

DO 3 I=1»N1 
TIMES! n = T( I) 

3 CONTINUE 
RETURN 

DETERMINE COEFFICIENTS FOR MODIFIED 
DEGREE ONE RATE FUNCTION 

4 TEST = -Al/ ( 2.=<‘A2 ) 

TINT = ER-EL 

IF (Al. GE. 0.. AND. A2.GT.0. ) GO TO 5 
GO TO 6 

5 B = A-A2=!‘TINT#=f'2 
Bl = Al+2.«A2*TINT 
GO TO 10 

6 B = A 

IF ( ( Al .LE. 0. .AND.A2.LT.0. ) . OR . ( Al . GT . 0 . . AN D . A2 . L T . 0 . 

AND. TEST. GE. TINT ) ) GO TO 7 
GO TO 8 

7 Bl = A1+A2*TINT 
GO TO 10 

8 IF (Al. GT. 0. . AND. A2.LT.0.. AND. TEST. LT. TINT) GO TO 9 
Bl = Al 

GO TO 10 

9 Bl = Al/2. 

MUST THE INTERVAL BE PARTITIONED? 

10 IF (Al*A2.LT.O. .AND. TEST. LT. TINT ) GO TO 11 
ERNEW = ER 

GO TO 12 

11 ERNEW = TEST+EL 

GENERATE DEGREE ONE NHPP ON INTERVAL 

12 BB = B 
BBl = Bl 

CALL NHPP2 (IS,EL,ERNEW,BB, BBl ,0,N1, lER) 

N( 1) = N1 

IF (N(l ).EQ.O) GO TO 14 

00 13 1=1, N1 
TIMES! I) = T( I ) 

13 CONTINUE 



COMPUTE LENGTH OF INTERVAL AND DETERMINE VALUE 
OF CSTAR for use in REJECTION ROUTINE 

14 Q = ERNEW-EL 
El = A 

E 2 = A2^Q*^2 
E3 = Al «Q 

E4 = Al=i==i'2/ (4.*A2 ) 

E5 = Al«=«'2/( 2.«A2) 

IF (Al .GE.0..AN0.A2.GT.0. ) GO TO 15 

IF (Al. LT.O.. AND. A2.GT.0.. AND. TEST. GE. TINT) GOTO 16 
IF (AI.LT.O.. AND. A2.GT.0. .AND.TEST.LT.TI NT) GO TO 17 
IF (AI.LE.0..AN0.A2.LT.0. ) GO TO 18 

IF (Al.GT.O.. AN0.A2 .lt. 0. .AND. TEST. GE.T INT) GO TO 19 
CSTAR = EXP(E1-E4)-EXP(E1) 

GO TO 20 

15 CSTAR = EXP(E1)-EXP(E1-E2) 

GO TO 20 
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16 CSTAR = EXP(E1+E2+E3)-EXP(E1+E3) 

GO TO 20 

17 CSTAR = EXP(E1-E4)-EXP(E1-E5) 

GO TO 20 

18 CSTAR = EXP( E1)-EXP(E1+E3+E2) 

GO TO 20 

19 CSTAR = EXP(EL+E3+E2)-EXP(E1) 

COMPUTE INTEGRAL OF MODIFIED DEGREE TWO RATE FUNCTION 
OVER INTERVAL 

20 CALL HELP ( A, A 1 , A2 ,EL , ERNE W t PMTR ) 

PMTR = PMTR-IEXP ( EXP( B1*ERNEW)-EXP(8 1=«'EL) ) ) /B1 

IDENTIFY AS FIRST SUBINTERVAL 
NOTE = 1 

GENERATE REALIZATION ON POISSON (PMTR) VARIATE 

21 CALL PVAR US, PMTR, M) 

IF (NOTE.EQ.l ) GO TO 22 
GO TO 25 

REJECTION ROUTINE USED ON FIRST SUBINTERVAL 

22 N(2) = M 
P(A) = B 
P(5) = 81 

IF (N(2).EQ.O) GO TO Ih 

CALL REJECT ( IS , EL ,CSTAR, P , Q, N ( 2) ) 

DO 23 I=1,M 
T IMES(N ( 1 ) + I) = T( I) 

23 CONTINUE 

HAS THE INTERVAL BEEN PARTITIONED? 

2A IF (ERNEW.EQ.ER) GO TO 34 
GO TO 27 

USE REJECTION ROUTINE ON SECOND PART OF INTERVAL 

25 N(4) = M 
P(4) = B 
P(5) = 81 

IF NO EVENTS OCCURRED BYPASS REJECTION ROUTINE 

IF (N(4),EQ.O) GO TO 35 
G = ER-ELNEW 

CALL REJECT ( I S , ELNE W ,C ST AR ,P , Q , N( 4) ) 

COPY TIMES OF EVENTS INTO ’TIMES* ARRAY 

N4 = N( 1 )+N (2 ) +N( 3 ) 

DO 26 I=1,M 
T IMES (N4+I ) = T( I) 

26 CONTINUE 

GENERATION OF VARIATES COMPLETE. 

GO TO ORDERING ROUTINE 

GO TO 35 

INTERVAL PARTITION WAS REQUIRED. MUST NOW 
CONSIDER SECOND SUBINTERVAL 

DETERMINE COEFFICIENTS FOR MODIFIED 
DEGREE ONE RATE FUNCTION 
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27 IF ( Al. GT.O. .AND. A2 .lt. 0. ) GO TO 28 
B = A-A2*TINT4'*2 

Bi = A1+2.«A2*TINT 
GO TO 29 

28 B = A+( AL/2 .)*TINT 
Bl = Al/2. + A2=:'TINT 

29 ELNEW = ERNEW 

GENERATE DEGREE ONE NHPP ON INTERVAL 

BB = B 
BBi = Bl 

CALL NHPP2 (IS,ELNEW,ERtBBtBBi»0,N3,lER) 

N(3) = N3 

IF (N(3). EQ.O) GO TO 31 
N3 = N( 1)+N(2) 

TRANSFER TIMES BETWEEN ARRAYS 

DC 30 1=1, N3 
TIMES(N3+I ) = T( I ) 

30 CONTINUE 

31 Q = TINT 

DETERMINE VALUE OF CSTAR FOR USE IN 
THE REJECTION ROUTINE 

E? = a^*Q*=!'7 
E3 = A1=4'Q 

IF (A1.GT.0..AND.A2.LT.0. ) GO TO 32 
CSTAR = EXP(E1-E4)-EXP( E1-E5-E3-E2) 

GO TO 33 

32 CSTAR = EXP(E1-E4)-EXP(E1+E3+E2 ) 

COMPUTE INTEGRAL OF MODIFIED DEGREE TWO RATE 
FUNCTION OVER SECOND INTERVAAL 

33 CALL HELP ( A, Al , A2 , ELNE W, ER ,PMTR ) 

PMTR = PMTR-(EXP(B)=«=(EXP(Bl*ER)-EXP(Bl*ELNEW) n/Bl 

IDENTIFY AS SECOND SUBINTERVAL 

NOTE = 2 
GO TO 21 

PARTITION OF INTERVAL NOT REQUIRED. COMPUTE TOTAL 
EVENTS AND SUPERPOSE TWO EVENT STREAMS 

34 N(5) = N(l)+N(2) 

IF (N(2).EQ.O) GO TO 38 
L8GN = N( 1 )+l 
JBGN = 1 

CALL COLATE ( L BGN ,N ( 5 ) , 1 ) 

GO TO 38 

PARTITION WAS REQUIRED. DETERMINE 
AMOUNT OF SORTING NEEDED 

35 N(5) = N( 1)+N(2)+N(3)+N(4) 

IF (N(2 ). EQ.O. AND. N(4 ). EQ.O) GO TO 38 
IF (N(4).EQ.O) GO TO 36 
IF (N(2).EQ.O) GO TO 37 

MUST SUPERPOSE FOUR EVENT STREAMS 

LBGN = N( 1 ) +1 
LFIN = N( 1)+N(2) 

CALL COLATE ( LBGN , LFIN, 1 ) 

LBGN = LFIN+N(3)+1 
JBGN = LFIN+1 

CALL COLATE ( LBGN ,N( 5 ) , JBGN ) 
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GO TO 3 8 

MLST SUPERPOSE FIRST HALF OF ARRAY ONLY 

36 N2 = N( 1)+N(2) 

LBGN = N( I )+l 

CALL COLATE (LBGN,N2,1) 

GO TO 38 

MUST SUPERPOSE SECOND HALF OF ARRAY ONLY 

37 KK = N( l)+N(2) + l 

LBGN = N(1 )+N (2)+N(3)+l 
LFIN = N(5) 

CALL COLATE ( LBGN ,N ( 5 ) » KK ) 

GO TO 38 

ARE TIMES OF EVENTS OR TIMES BETWEEN EVENTS REQUESTED? 

38 IF (lI.EQ.O) RETURN 

CALCULATE TIMES BETWEEN EVENTS 

S = TIMES( 1 ) 

TIMES! 1 ) = TIMES! l)-EL 
N5 = N( 5) 

DO 39 1=2, N5 
SI = TIMES! n 
TIMES!! ) = TIMES! I )-S 
S = SI 

39 CONTINUE 
RETURN 
END 



SUBROUTINE NHPP2 SIMULATES A NON-HOMOGEN EOU S 
PCISSON PROCESS WITH A LOG-LINEAR INTENSITY 
!RATE) FUNCTION 

SUBROUTINE NHPP2 ( I S , EL , ER , A, A 1 , 1 1 , N , I £R ) 

DIMENSION T!5000) 

COMMON /HOLD/ T 

CALL OVFLOW 

INITIALIZE VARIABLES 

lER = 0 
TINT = ER-EL 
A = EXP!A+A1*EL) 

IS THE POISSON PROCESS HOMOGENEOUS? 

IF !A1.EQ.0.) GO TO 3 

PAR = ! A’l' ! EXP!TINT*A1 )-l. ) )/AI 

IF !A1.GT.0.) GO TO I 

IFLAG = 3 

GO TO 2 

1 A = A=<=EXP! TINT*Al ) 

Ai = -A1 

IFLAG = 2 

COMPUTE parameters OF BOTH POISSON RANDOM VARIABLES 

2 THETA = -A/Al 
GC TO 4 

COMPUTE RATE AND SCALING PARAMETERS FOR HOMOGENEOUS 
PCISSON PROCESS 
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3 PAR = TINT^i^A 
IFLAG = I 
AINVRS = l./A 

COMPUTE NUMBER OF EXPONENTIAL VARIATES REQUIRED 

4 NMAX = PAR+6.*SQRT(PAR) 

IS THIS A HOMOGENEOUS POISSON PROCESS? 

IF (IFLAG. EQ.i) GO TO 17 

GENERATE REALIZATION ON POISSON (THETA) VARIATE 

5 CONTINUE 

CALL PVAR (IS, THETA, M) 

IF (M.EQ.O) GO TO 7 

CALCULATE TIMES OF EVENTS 

CALL SEXPON (IS, T, NMAX) 

B = -A1 

V = 0. 

UMAX = NMAX+1 
DC 6 1 = 1, UMAX 

HAVE NUMBER OF EVENTS EXCEEDED THE MAXIMUM NUMBER 
THAT THE ARRAY CAN HOLD? 

IF ( I.GT.NMAX) GO TO 8 

V = V+T(I ) /((M-I+l )*B) 

IF (V.GT.TINT) GO TO 9 
T ( I ) = V 

IF (I.EQ.M) GO TO 10 

6 CONTINUE 

NO EVENTS OCCURRED 

7 N = 0 
RETURN 

TOO MANY EVENTS FOR ARRAY. INCREMENT ERROR 
CODE AND TRY AGAIN 

8 lER = lER+l 
GO TO 5 

THE NUMBER OF EVENTS OBSERVED TO OCCUR IN THIS 
NON-HOMOGENEOUS POISSON PROCESS IS ‘N’ 

9 N = I-l 

IF (N.EQ.O) RETURN 
GO TO 11 

10 N = M 

11 CONTINUE 

IS THE RATE FUNCTION INCREASING OR DECREASING? 

IF (IFLAG. EQ.3) GO TO 13 

TIME REVERSAL TECHNIQUE IS NECESSARY 
DETERMINE WHETHER N IS EVEN OR ODD 

I SIG = M00( N, 2 ) 

NLOOP = N/2 
N1 = N+1 

00 12 1=1, NLOOP 
S = T ( I ) 

T( I) = ER-T(N1-I ) 

T(Nl-I) = ER-S 
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12 CONTINUE 

IF (ISIG.EQ.l) T(NL00P+1)=ER-T(NL00P+1) 

ARE TIMES OF EVENTS REQUESTED? 

IF (lI.EQ.O) RETURN 
GO TO 15 

13 IF ( II.NE.O) GO TO 15 
IF (EL.EQ.O.) RETURN 

DO 14 1=1, N 
T( I) = EL+T(I ) 

14 CONTINUE 
RETURN 

CALCULATE TIMES BETWEEN EVENTS 

15 S = T(1 ) 

DO 16 1=2, N 
SI = T ( I ) 

T( I) = T( I )-S 
S = SI 

16 CONTINUE 
RETURN 

THE POISSON PROCESS IS HOMOGENEOUS 

17 I = 1 
U = 0. 

CALL SEXPON (IS,T,NMAX) 

18 0 = U+T( I ) 

IF (U-GT.PAR) 30 TO 20 
I = I+l 

IF (I.GT.NMAX) GO TO 19 
GO TO 18 

INCREMENT ERROR CODE 

19 lER = I ER+1 

TRY AGAIN WITH NEW STRING OF VARIATES 

GO TO 17 

20 N = I-l 

IF (N.EQ.O) RETURN 
IF (lI.EQ.l ) GO TO 22 

DO 21 1=1, N 

EL = EL+AlNVRS^i'T ( I ) 

T( I ) = EL 

21 CONTINUE 

RETURN 

22 DO 23 1=1, N 

T( I) = T( I )*AINVRS 

23 CONTINUE 
RETURN 
END 



SUBROUTINE PVAR GENERATES A POISSON (THETA) 
VARIATE, M, USING THE GAMMA METHOD 

SUBROUTINE PVAR (IS, THETA, M) 

DIMENSION T(5000) 

COMMON /HOLD/ T 
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K = 0 
C = 16. 

D = . 875 

1 IF (THETA. LT.C) GO TO 2 
GO TO 5 

2 U = 1. 

CTN = EXP(-THETA) 

MMAX = THETA+6.*SQRT(THETA) 

3 1 = 1 

CALL SRANO (IS»T, MMAX) 

4 U = U*T(I) 

IF (U.LT.CTN) GO TO 8 
I = I + l 
K = K+1 

IF (I.GT.MMAX) GO TO 3 
GO TO 4 

5 NP = INT(D*THETA) 

AN = FLOAT(NP) 

CALL GAMA ( ANt IS» G, 1 ) 

IF (G.GT. THETA) GO TO 6 
K = K+NP 
THETA = THETA-G 
GO TO 1 

6 U = THETA/G 
NP = NP-1 

CALL SRANO (IS»TtNP) 

DO 7 I=l»NP 

IF (T(I ) .LT.U ) K = K+1 

7 CONTINUE 



THE VALUE M IS ASSUMED BY THE POISSON (THETA) VARIATE 

8 M = K 
RETURN 
END 



SUBROUTINE REJECT GENERATES AN ORDERED SAMPLE 
OF GIVEN SIZE FROM THE SECOND COMPONENT 
OF THE ORIGINAL INTENSITY FUNCTION 
USING A REJECTION-ACCEPTANCE ALGORITHM 

SUBROUTINE REJECT ( I S , EL ,CS TAR , P VEC , Q t L ) 
DIMENSION V(500), PVEC(5) 

DIMENSION T(5000) 

CCMMON /HOLD/ T 
L20 = L*10 

IF (L20.GT.500) L20=500 
LI = L+1 
K = 1 

1 J = 0 

CALL SRANO (IS,V,L20) 

DO 2 I=1tL20 
J = J + 1 

T(K) = Q#V(J)+EL 
J = J+1 

IF ( V( J ) .LT .CALC ( PVECtT (K) )/CST AR) K=K+1 
IF (K.EQ.Ll) GO TO 3 
IF ( J.GE.L20-1 ) GO TO 1 

2 CONTINUE 

IF (K.LT.L) GO TO 1 

3 CALL PXSORT (T, 1 » L) 

RETURN 

END 
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SUBROUTINE HELP EVALUATES THE INTEGRATED 
FUNCTION OVER THE INTERVAL (EL,ER) 

SUBROUTINE HELP { A , A1 , A2, EL » ER , XX ) 

DOUBLE PRECISION MMOAW,BBtAA 
Z = SQRT( ABS(A2) ) 

Y = (AL^S'Z)/ (2.*A2) 

AA = Z=s=EL+Y 
SB = Z=<'ER+Y 
CC = A-Al^i'Al/I^.^AZ) 

CC = EXP(CC)/Z 

IF (A2.LT.0.) GO TO 1 

Q1 = DEXP( AA=i'=<'2 ) =*'MMDAW( AA ) 

Q2 = 0EXP( BB^^i'ZI + MMOAWIBB) 

XX = CC*(Q2-QL) 

RETURN 

CC = CC*. 8862269 

XX = CC*(OERF(BB)-DERF(AA ) ) 

RETURN 

END 



INTENS ITY 



SUBROUTINE COLATE SUPERPOSES TWO ORDERED 
EVENT STREAMS OVER THE SAME INTERVAL 

SUBROUTINE COLATE { LBGN , LF I N , J BGN ) 
DIMENSION TIMES(5000), T(5000) 

COMMON /MIKE/ TIMES/HOLD/T 
I = JBGN 
J = I 
K = LBGN 

1 IF (TIMES! I ).LT.TIMES(K) ) GO TO 2 
T(J) = TIMES(KJ 

J = J+1 
K = K + 1 

IF (K.GT.LFIN) GO TO 3 
GO TO L 

2 T(J) = TIMES(I) 

J = J+ 1 

I = I + l 

IF ( I.EQ.LBGN) GO TO 5 
GO TO L 
311= LBGN-1 

DO 4 N= 1 , 1 1 
T(J) = TIMES(N) 

J = J+l 

4 CONTINUE 
RETURN 

5 CONTINUE 

DC 6 N=K,LFIN 
T(N) = TIMES(N) 

6 CONTINUE 

RETURN 

END 



FUNCTION CALC EVALUATES THE SECOND COMPONENT OF 
THE DECOMPOSED INTENSITY FUNCTION 
FOR ANY INPUT VALUE. 

FUNCTION CALC (P,ABSA) 

DIMENSION P(5) 

X = P( 1 )+P (2)*ABS A+P(3 >*A3SA**2 
XX = P( 4)+P(5)*ABSA 
CALC = EXP( X)-EXP (XX) 

RETURN 

END 
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